Lattice imperfections and high-harmonic generation in correlated systems

We study effects of lattice imperfections on high-harmonic generation from correlated systems using the Fermi–Hubbard model. We simulate such imperfections by randomly modifying the chemical potential across the individual lattice sites. We control the degree of electron–electron interaction by varying the Hubbard U. In the limit of vanishing U, this approach results in Anderson localization. For nonvanishing U, we rationalize the spectral observations in terms of qualitative k-space and real-space pictures. When the interaction and imperfection terms are of comparable magnitude, they may balance each other out, causing Bloch-like transitions. If the terms differ significantly, each electron transition requires a relatively large amount of energy and the current is reduced. We find that imperfections result in increased high-harmonic gain. The spectral gain is mainly in high harmonic orders for low U and low orders for high U.


Introduction
Studying the dynamics of solids on the ultrafast (10 −18 s-10 −15 s) timescale is an interesting possibility.One of the approaches used to explore these timescales is high-order harmonic generation (HHG) [1][2][3][4][5][6].In addition to exploring dynamics, HHG is used to produce ultrashort pulses with high photon energy [1,[7][8][9][10].The high density of solids may aid in the production of intense HHG pulses and thereby improve their usefulness for time-resolved measurements [5,7,[11][12][13][14][15][16][17][18][19][20][21].Describing HHG from bandgap materials is often done within the framework of bandstructure theory.Here HHG has been qualitatively separated into two categories, intra-and interband HHG.Intraband HHG is the generation of harmonics by the acceleration of charge carriers within a given band.Interband generation is perhaps best described using the so-called 3-step model: (1) an electron is excited across the bandgap, (2) the electron and hole propagate through their respective bands, resulting in intraband HHG, (3) the electron and hole recombine under the emission of an interband photon with energy corresponding to the difference between the bands at the instant of recombination [22,23].The second step includes charges propagating through bands, resulting in intraband HHG.The two mechanisms are therefore intrinsically linked.Furthermore, the total HHG spectra arise from a coherent sum of the two contributions to HHG which is why the separation is merely qualitative [24,25].
The band picture used for the 3-step model requires the translational symmetry gained by applying the mean-field approximation.The mean-field approximation allows the electrons to be treated as independent particles, which is not exact.With that in mind, over the last 5 years, researchers have worked on the effects of beyond mean-field electron-electron interaction on HHG both theoretically [6,21,[26][27][28][29][30][31][32][33][34][35][36][37] and experimentally [38][39][40].Much theoretical research [6,21,27,33,37,41] has used the Hubbard model, which is also used in this work.This model includes electron-electron interactions through the Hubbard U, which energetically penalizes the occurrence of doubly-occupied sites.Studies using the Hubbard model have led to the development of a parallel to the 3-step model but in the quasiparticle picture of doublons and holons [21].Doublons are doubly occupied lattice sites, and holons are empty lattice sites.Mott insulators, i.e. solids dominated by electron-electron interactions and with an uneven number of valence electrons per lattice site, are typically well described by this quasiparticle picture.The doublon-holon 3-step model consists of the following steps, (1) a doublon-holon pair is created, (2) the doublon and holon propagate through the solid, and (3) they recombine [21,26,37,42].The onsite interaction term, U, causes the creation of a doublon to require a significant energy, equivalent to the excitation of an electron across a band gap.Similarly, the recombination of a doublon and holon corresponds to the recombination of the electron and hole from the original 3-step model.Other studies of HHG from correlated systems have led to related descriptions.An example is a study of bipartite lattices where the doublon-holon excitations are replaced with a pseudo-spin kink and anti-kink type excitation in a 3-step-type model [29].In the present work, we include only onsite U, but inter-site electron-electron repulsion can affect the doublon-holon picture described above, see, e.g.[32].
The Hubbard model describes a translationally symmetric lattice, which is an idealization as lattice vibrations cause aperiodicity in the lattice.Moreover, defects and impurities break the translational symmetry and have visible effects on HHG spectra, see, e.g.[43] and references therein.Likewise the presence of donors locally changes the potential energy and destroys the lattice translational symmetry and these effects affect the harmonic spectra too [44].To capture some essential elements of such lattice irregularities in the simplest possible setting while keeping the interactions between the electrons, we introduce an extra term to the Hubbard Hamiltonian.This term seeks to encapsulate the above physics by modelling the essential randomness by requiring the energy to vary from site to site.It was shown long ago in the seminal work of Anderson [45] that such alterations of the potential energy landscape significantly influences the system conductivity by leading to the localization of the wave functions at the sites with most deep potentials, and hence reducing the conductivity, a phenomena now known as Anderson localization.In connection with HHG, studies have been performed from lattices with imperfections for U = 0 [43,[46][47][48][49][50][51][52][53].
In those studies, Anderson localization caused high-harmonic gain [51].The study of impurities has led to a 3-step model qualitatively similar to that of atomic HHG.In the first step, an electron or hole is excited away from the impurity, in the second step, the unbound electron or hole propagates within the lattice, and in the third step, it recombines with the impurity resulting in the emission of high harmonics [46].
Our present work addresses fundamental scientific issues and elucidates the following research questions for the first time: (i) How are HHG spectra affected by having a lattice with a degree of randomness in the site-dependent potentials in the presence of electron-electron interactions?(ii) How do effects of having such randomness in the lattice and electron-electron interactions relate?and (iii) How do dynamics in the lattice relate to the spectra?Even though experiments on correlated materials subject to random insertion of donors or lattice distortions are not available at the time being, they may appear in the future, and our work would be relevant for interpretation and physical insight.
The paper is organized as follows.In section 2, we introduce the model and methods.In section 3, we discuss two qualitative physical pictures that are then used to discuss our results.In section 4, we summarize and conclude.The appendix explores the sensitivity of the results to the sampling method used for the simulation of the imperfections.Atomic units are used throughout unless stated otherwise.

Hamiltonian
We use the Hubbard model to simulate a one-dimensional L = 12 site chain with periodic boundary conditions at half filling, i.e. the number of electrons equals L. Representing a solid in this manner implies assuming that the dynamics at one site is unaffected by the dynamics at sites more than L sites away.Therefore we may still consider periodic boundary conditions when we add a term to model lattice imperfections.Treating the laser pulse in the dipole approximation results in the Hamiltonian with the laser-driven Hubbard Hamiltonian given by where In equation (1), Here ĤHop (t) describes hopping, ĤU describes electron-electron interactions, and Ĥ∆ treats the lattice imperfections.Within ĤHop (t), −h 0 is the transition matrix element corresponding to nearest-neighbor hopping, e iaA(t) is Peierls phase describing the laser-matter interaction, and c † i,σ and c i,σ are creation and annihilation operators of an electron on site i with spin σ ∈ {↑, ↓}.In Peierls phase, a is the lattice spacing and A(t) is the vector potential of the laser pulse.In ĤU , U describes beyond-mean-field energy effects resulting from a doubly occupied site, i.e. a doublon, and ni,σ = ĉ † i,σ ĉi,σ counts the number of electrons on site i with spin σ.Finally, in Ĥ∆ , ∆A i is the energy needed to add an electron to site i.We pick the A i values according to equation ( 6) by evaluating cosine of a random, uniformly distributed phase ϕ i ∈ [0; 2π], and finally in equation ( 5) scale A i by the system parameter ∆.It is common to generate the A i -values with this distribution, see, e.g.[51,52].As the evaluation of the time-dependent Hubbard model is numerically challenging even without introducing extra effects, such as disorder, we choose the above modelling of the disorder.This methodology was, e.g.previously used to study quasicrystals [52], which show similar lattice structure as one would expect from a system with many phonon modes activated [54].The energy axis is set such that ⟨A i ⟩ = 0.All simulations use the same ϕ i -values unless we explicitly state otherwise.Note that picking A i with alternating signs on the different sites yields the ionic Hubbard model [55,56].For nonvanishing ∆, the Hamiltonian in equation ( 1) shifts the energy on sites in the lattice that are occupied by electrons, as seen from the expression on the right-hand side of equation (5).Clearly these site-dependent shifts destroy the translational symmetry of the system.In a first-quantization formalism the site-dependent shift of potential energy means that the local potential, say V(x), varies from site to site, and then V(x + a) is of course not equal to V(x), and the periodicity is lost.
Ordinarily, in perfect lattices, the onsite matrix element is set to zero by adjusting the zero-point energy accordingly [6,21].This choice effectively sets the chemical potential to zero.Ĥ∆ 's effects can be described, in physical terms, as modifying the chemical potential at each site by ∆A i .It is important to stress that the present model only considers changes in the local potential energy-if an electron occupies a site, it is ascribed a random energy, as described by equations ( 5) and (6).When we use the wording 'lattice imperfections' in the remainder of the paper, it is this kind of random distortion that we consider, and no changes of the lattice constant a are considered.For a review on localization phenomona including a discussion of lattice models, see, e.g.[57].

Measures
We consider a system that is thin in the direction of the laser propagation.This thinness allows us to neglect propagation effects [58].The spectrum, S(ω), is then proportional to the norm-square of the Fourier transform of the electron acceleration [59] where j(ω) is the Fourier transform, F, of the current j(t).The current can be determined via Heisenberg's equation of motion [60], resulting in

Simulation method and parameters
We solve the time-dependent Schrödinger equation i d dt |Ψ(t)⟩ = Ĥ(t)|Ψ(t)⟩ for |Ψ(t)⟩ expressed in terms of configurations |Φ i ⟩, with i being a collective index specifying the distribution of all the electrons on the sites, i.e. |Ψ(t)⟩ = i c i (t)|Φ i ⟩.The simulations use the Arnoldi-Lanczos algorithm with Krylov subspace decomposition, see, e.g.[61][62][63][64].The ground state of the system is chosen as the initial state in our simulations.We obtain the ground state by propagation in imaginary time, using a time-step size of 1/ √ 10 and a maximal Krylov subspace dimension of 4, which will also be used for the real-time propagation.The imaginary-time propagation is taken to be converged when the change in energy of the state over 100 time steps is limited by machine precision.The real-time propagation convergence is tested by comparing spectra and currents for different Krylov subspace dimensions and time steps.We choose the system parameters to describe Sr 2 CuO 3 [65], as was done previously [6].The specific values are a = 7.5589 a.u. and h 0 = 0.0191 a.u.The laser pulse is an N c = 10 cycle pulse with a vector potential given by Here, the laser angular frequency is ω L = 0.005 a.u. and the vector potential amplitude A 0 = F 0 /ω L = 0.194 a.u.corresponds to a peak intensity of 3.3 × 10 10 W cm −2 .The Hubbard U and the disorder strength ∆ are treated as parameters.

results
In this section, we start by discussing two qualitative pictures, which can be used as frameworks for rationalizing the results.

Qualitative physical pictures 3.1.1. k-space quasiparticle picture: interaction based dynamics
We employ the doublon-holon quasiparticle picture as a framework to understand the HHG spectra from the highly correlated perfect lattice.In the highly correlated limit the creation of doublons requires a significant amount of energy.This observation means that the energy of a given eigenstate is closely connected to the number of doublons one would expect to observe in that eigenstate.Note that as h 0 is non-zero the eigenstates generally contain configurations with a varying number of doublons, but with the majority of the population in configurations with a specific number of doublons.This means the eigenstates are grouped into states with similar energy and doublon count.These groups are a reflection of the Hubbard sub-bands, also known as doublon-holon bands [21,37].In figure 1 we show a histogram of the eigenenergies of the U = 10h 0 , ∆ = 0 system, which were found by numerical diagonalization.Note that the correlated system without lattice imperfections can be diagonalized using the Bethe ansatz [66].From this approach a quasiparticle bandstructure can be derived.However, Ĥ∆ breaks the symmetry of the system and the quasiparticle bandstructure picture breaks down for ∆ > 0. The histogram shows seven distinct groups of states, each associated with a given number of doublons.Notice the two small groups, (i) below the zero energy and (ii) between 230 and 240 ω L of energy.In a lattice with L = 12 sites and L electrons, there are electron configurations containing between zero and six doublons, corresponding to seven groups of states.
The band structure generally consists of bands with bandwidth very close to the width of the Bloch band, described by ĤHop (t), and separated by a bandgap, induced by ĤU .This band gap results from the energy penalty associated with the creation of doublons, i.e.U.
As the system is half-filled, having zero doublons in an electron configuration requires having precisely one electron at each site.This observation means that any single-electron transition from such a configuration will result in a doublon-holon pair.As configurations with one electron on each site are dominant in the ground state of the system, it is exceedingly unlikely for the system to transition from the ground state to another state in the lowest group.This conclusion means that a rather significant amount of energy is required to drive the system out of the ground state.That amount of energy is commonly referred to as the Mott-gap and is indicated by the blue arrow in figure 1.The Mott-gap is commonly calculated as, Here E N GS (U, ∆) is the ground state energy of the system with N electrons.Note that half-filling corresponds to N = L.We highlight the ground state energy dependence on the U-and ∆-variable explicitly.This equation can be justified as follows.Adding one electron to an otherwise half-filled system results in an increase in energy.This energy increase comes mainly from the doublon created by adding an extra electron.The remainder of the energy comes from the Bloch-band described by ĤHop (t).Removing an electron from the half-filled system results in the removal of a bit of energy, equivalent to the energy from the Bloch-band added by the addition of an electron.This means E L+1 GS (U, ∆) + E L−1 GS (U, ∆) is essentially equal to the energy required to create a doublon plus twice the ground state energy of the half-filled system.So by subtracting twice the ground state energy of the half-filled system, ∆ Mott (U, ∆) becomes a decent approximation of the amount of energy required to create a doublon, i.e. the bandgap between the Hubbard-sub-bands.
Once the system has been excited to include a doublon-holon pair those quasiparticles may propagate through the lattice.In doing so, the system propagates through the Hubbard-sub-band in question resulting in intra-Hubbard-sub-band HHG.At some later point, they may recombine resulting in inter-Hubbard-sub-band HHG.This constitutes another way to view the doublon-holon 3-step model  10), and the green arrows indicate the width of the Bloch band, ∆ width of equation (11).mentioned in section 1.Firstly, the system is excited across the bandgap between the Hubbard sub-bands, under the creation of a doublon-holon pair.Secondly, the system propagates through the sub-band, as the doublon and holon propagate through the lattice and, finally, the system drops back to the lower-lying bands, through recombination of the doublon-holon pair.The minimum energy released from such a recombination event is the Mott gap, and the maximum energy is the Mott gap plus twice the width of the Hubbard sub-band, as the state may fall from the top of the upper band to the bottom of the other.
As the width of the Hubbard-sub-bands is very similar to the width of the Bloch bands, we have indicated the Bloch-band width with the green arrows in figure 1.The bandwidth of the Bloch band is which follows from the tight-binding Hamiltonian ĤHop (t).We, therefore, expect the inter-Hubbard-subband HHG to dominate between the Mott-gap, i.e.
and the Mott gap plus two Bloch-band widths, i.e.
The intra-Hubbard-sub-band HHG can at most result in harmonics with energy corresponding to the Hubbard-sub-band width, i.e. essentially the Bloch bandwidth.We will use these energy limits in the interpretation of the HHG spectra in section 3.

Real-space picture: imperfection-induced dynamics
We now turn to a qualitative discussion of the effects of imperfections.In the limit of h 0 ≪ max(U, ∆), the eigenstates of the systems converge toward individual electron configurations.As those configurations are the basis states used in the simulations, finding the ground state is trivial in this limit.Analysis of those ground states gives insight into the possible dynamics of the system.We have illustrated such dominant ground state configurations in figure 2 for four different values of ∆ in terms of U for the specific realization of the set of {A i } values used in the majority of the simulations.Whether one considers the change of ∆ in reference to U or views it as a change in U, the figure, results, and conclusions will be identical.Notice the significantly varying ordinate scale in the panels.In the limit of U ≫ (∆, h 0 ), the system is a Mott insulator due to the half-filling of the system.For ∆ = 0.1U, figure 2(a), the effect of Ĥ∆ is almost insignificant, notice the scale of ordinate.As a result, ĤU dominates the system, which results in one electron on each site and arranged in antiferromagnetic ordering.This ordering means any transition across nearest-neighbor sites creates a doublon-holon pair, which requires an energy ∼U.For sufficiently large U, compared to ∆, that energy requirement can effectively lock all the electrons resulting in ineffective harmonic generation [6,21,37,41].In this limit, the doublon-holon-based 3-step model captures the physics as ∆ has minimal effect.For ∆ = U, figure 2(b), Ĥ∆ is significant enough to create three doublons, on sites 4, 6, and 7.It is now, e.g.significantly easier for the electron on site 10 to hop to either site 9 or 11 than in the ∆ = 0.1U case.As ∆ increases, the doublon count in the ground state increases.We see this trend in figure 2(c), where ∆ = 4U, and sites 0, 4, 6, 7, and 10 contain doublons, and finally in figure 2(d), where ∆ = 10U, and all the electrons are in doublons on sites, 0, 2, 4, 6, 7, and 10.This trend illustrates the interplay between Ĥ∆ and ĤU , where Ĥ∆ tends to collect the electrons on the sites with the lowest ∆A i -values, and ĤU tends to push the electrons apart, resulting in an even electron distribution.Describing this in terms of the impurity-based 3-step model, the increasing ∆ corresponds to the lowered binding energy, which could increase the maximal harmonic orders created but at the same time lower the gain from the given imperfection.
It is also worth paying attention to the transition from site 7-8.In figure 2(d), ∆ = 10U, exciting an electron from site 7-8 would require ≃ 18.2U = 1.82∆ of energy, whereas for the smaller ∆ = 4U, seen in figure 2(c), it only requires 6.68U = 1.67∆ and for ∆ = U, seen in figure 2(b), the gap is down to 0.92U = 0.92∆, finally for ∆ = 0.1U requires 1.19U = 11.9∆.Therefore, lowering ∆ in terms of U could decrease the energy required for specific transitions.When ∆ ≫ U the binding energy is massive, it then falls as U increases until the U term results in a more even electron distribution across the lattice, meaning the doublon-holon-based 3-step model becomes the appropriate picture.Finally, we point out the dynamics as those associated with site 5 in figures 2(b) and (c).The addition of doublons to the ground state increases the amount of Pauli blocking in the system, which causes the electron on site 5 to be frozen in place, due to the doublons on sites 4 and 6 and the nearest-neighbor hopping approximation.We note that this effect operates on a site-to-site basis, which is not experimentally testable, and appears to have little impact on the HHG spectra.Therefore, we will not comment further on this effect in this work.We also note that in higher dimensions the number of sites with comparatively low A i -values would rise significantly making this effect quite unlikely, which further decreases the importance of this effect.In conclusion, it is possible, as illustrated in figure 2, for the electron-electron repulsion of ĤU to partially balance out the attraction to specific sites induced by Ĥ∆ .That balancing can make it easier to drive electrons across sites.This mechanism works even for large U, when ∆ is large enough.Note also that this balancing act makes it virtually impossible to tell intra and inter doublon-holon band HHG apart.Telling the two types apart requires that doublon-holon bands are well defined, which is not the case when ∆ ≃ U. Therefore, in this limit, any observed enhancement resulting from a significant ∆ value cannot be directly attributed to either mechanism.

∆ Mott (U, ∆)
Before presenting the HHG spectra, we consider the dependence of the Mott gap on the magnitude of ∆ and the choice of A i -values, as this dependence supplements the physical picture discussed in section 3.1.2.In figure 3 we show ∆ Mott (U = 10h 0 , ∆ = nh 0 ) for all integer values of n between 0 and 10, and 10 different sets of A i -values.Each set of A i -values corresponds to a given color.We only calculate ∆ Mott (U, ∆) for U = 10h 0 as the Mott gap is only descriptive in the high-U limit, and U = 10h 0 is the highest U used in this work.The figure shows that the Mott gap decreases with ∆ up to ∆ ≃ 5h 0 , and then the gap stays at a constant level of around 0.015 with increasing fluctuations for different {A i } sets.Firstly, ∆ Mott (U = 10h 0 , ∆ = 0) does not depend on the choice of A i -values, see equation (5).We observe a change in the slope at ∆ ≃ h 0 .When ∆ ≲ h 0 , ∆ Mott (U = 10h 0 , ∆) appears to drop relatively slowly.This can be understood through the qualitative picture of section 3.1.2,as well as the text surrounding equation ( 10): When U is large the ground state of the half-filled system essentially contains 0 doublons.This means adding an electron creates a doublon in the ground state, see discussion of equation (10).This doublon is distributed evenly across the lattice when ∆ = 0, due to the ĤHop (t) term.This even distribution means that for low ∆ the energy of the doublon is less affected by changes in ∆ as ⟨A i ⟩ = 0. Similarly, the hole created by removing an electron from half-filling starts evenly distributed, which also means the associated energy is largely independent of ∆.This means that for ∆ ≈ 0, ∆ Mott (U, ∆) is less affected by changes in ∆, than otherwise for ∆ < 5h 0 , see figure 3 for the lowest ∆ values.As ∆ increases the doublon and hole localize to the sites with minimal and maximal A i values, respectively.This localization causes the scaling of ∆ Mott (U, ∆) to become approximately linear in ∆ for ∆ between h 0 and 5h 0 .As the effects of ĤHop (t) on ∆ Mott (U, ∆) becomes largely insignificant, the ∆ scaling is determined from Ĥ∆ , which scales linearly in ∆.
When ∆ goes beyond 5h 0 , ∆ Mott (U, ∆) fluctuates around 0.01-0.02.Note that as the A i values are distributed between -1 and 1, Ĥ∆ can create an energy difference between sites of, at most, 2∆.This means that the energy difference in Ĥ∆ from removing or adding an electron is at most 2∆, whereas from ĤU it is merely U. Thus the maximal induced energy difference from Ĥ∆ equals that from ĤU when ∆ = U/2 = 5h 0 .So at this point the effects of ĤU and Ĥ∆ become comparable and the doublon-holon picture of the dynamics described in section 3.1.1breaks down.When ∆ is so large, the ground state of the half-filled system is likely to contain doublons, as explained in section 3.1.2.As a result of these doublons the reasoning behind the equation for ∆ Mott (U, ∆), equation ( 10), ceases to be accurate, as one can no longer claim that adding an electron to the lattice will create a doublon.This means that ∆ Mott (U, ∆) largely loses its relevance.We do note that the variance in ∆ Mott (U, ∆) over sets of A i values seems to scale approximately linearly with ∆ in this regime, which is to be expected given the form of Ĥ∆ .

Spectra
We now discuss the main results of this work, i.e. the high-harmonic spectra.In figure 4, HHG spectra for various U and ∆ values are shown.In all panels, the energy width of the Bloch band (equation ( 11)) in units of the laser angular frequency is shown by full vertical lines at ∆ width /ω L ≃ 15.3.The leftmost vertical dashed lines in panels (b), (c) and (d), indicate the Mott gap of equation ( 12), while the rightmost vertical dashed lines indicate the value of the Mott gap including twice the Bloch-band bandwidth of equation (13).In previous studies with ∆ = 0, it was shown that inter doublon-holon band HHG dominates between the dashed vertical lines in figure 4(d), i.e. between the Mott gap and the Mott gap plus twice the Bloch band width.In the Mott limit, intra doublon-holon HHG tends to be the main source of HHG below the Mott gap [21,37].In figure 4(a) [U = 0] the ∆ = 0, Bloch intraband result is shown in dark blue.As Ĥ∆ consists of one-particle operators describing electron-imperfection scattering, a small ∆ only facilitates harmonics up to the bandwidth of the Bloch band corresponding to about fifteen harmonic orders.We note that the ∆ = 1h 0 , ∆ = 2h 0 , and ∆ = 5h 0 all begin to drop-off around that harmonic order.In the ∆ = 10h 0 case, the spectrum shows several significant localized enhancement features, e.g. at harmonic orders around 39 and 58.These arise from specific energy difference induced in the system by Ĥ∆ , and we will return to this point below.
As U increases relative to ∆, the electrons become bound in antiferromagnetic ordering and the doublon-holon picture (section 3.1.1)of HHG becomes accurate, see also [21,37].In short ∆ Mott (U, ∆) acts as a bandgap, meaning a parallel to interband HHG is observed at harmonic orders between ∆ Mott (U, ∆)/ω L (equation ( 12)) and (∆ Mott (U, ∆) + 2∆ width )/ω L (equation ( 13)).This results in a plateau-like structure between the two values, i.e. between the dashed black lines in figure 4.This plateau is generated mainly from inter doublon-holon band transitions.If ∆ Mott (U, ∆) > ∆ width , then intra doublon-holon band HHG, that is HHG from doublon-holon pair preserving transitions, is also visible in the spectra, mainly below the 15 harmonic orders corresponding to ∆ width /ω L .Within this plateau the spectra are rather irregular.This irregularity been seen in previous works using similar systems, e.g,.[6,21,26,33,37,67], and has been speculated to originate from limited pulse lengths and missing dephasing channels [21].In [6], the less clear peaks at the odd harmonics were attributed to the U term.If both ∆ and U are large, then it becomes possible to have a doublon situated at an imperfection in the ground state, see figure 2. Ordinarily, large ∆ would imply large binding energy and, therefore, high harmonic orders.However, large U reduces that binding energy significantly resulting in high harmonic gain for low harmonic orders with increasing ∆ as observed in figure 4(d).Note that in the ground state of the high-U and low-∆ systems, any dynamics in the system must begin with the creation of a doublon-holon pair, i.e. an inter doublon-holon band transition.For the intermediary U values in figures 4(b), U = 2h 0 , and figure 4(c), U = 5h 0 , the U value is not so large as to freeze the dynamics significantly, but large enough that it impacts the binding of electrons to the lattice imperfections.In figure 4(b), this balance manifests itself as very little change in the spectra regardless of the degree of impurity.This lack of change likely results from the amount of energy required to create a doublon being similar to the energy required to remove an electron from a binding imperfection.In figure 4(c), the differences between the spectra with different ∆ values are more pronounced than in figure 4(b).In figure 4(c), U affects the binding energies significantly.These binding energies mean the energies involved in HHG from impurities are lowered, compared to figure 4(b), while the higher U means more energy is required to produce a doublon.This observation explains the trend towards higher gain for lower ∆ values.The overall lower gain in the high-∆ cases, compared to the ∆ = 0 result, is likely a result of the fewer sites that facilitate the impurity-based HHG compared to the doublon-holon-based HHG.HHG from impurities can only originate from sites with binding impurities whereas doublon-holon pair creation is site independent.Another way to describe the change in these spectra, figures 4(b) and (c), is by considering the effect of Ĥ∆ on the doublon-holon picture of section 3.1.1.The ∆ term effectively modifies the energies associated with the interband transitions.This means Ĥ∆ effectively softens the plateau-like structure of the inter-doublon-holon-band HHG.This is clearly seen in figure 4(c) as the observed drop in intensity with increasing ∆.In figure 4(d), U = 10h 0 , this softening effect is not visible due to the significant enhancement for low harmonic orders.This enhancement arises due to the ∆-induced lowering of the energy required to generate a doublon-holon pair, and is therefore not necessarily an enhancement of intra doublon-holon band HHG despite the low harmonic orders.That enhancement leads to a drop-off that does not reach the intensities of the inter-doublon-holon-band HHG before around the fifty-fifth harmonic, i.e. around the harmonic corresponding to (∆ Mott (U, ∆) + 2∆ width )/ω L .

Effects of specific A i -values
We now turn to a discussion of the localized enhancement features mentioned earlier.These enhancement features are related to the specific set of disorder variables used, i.e. the specific set of ∆A i values in equation (5).As an example, we consider in figure 5 the spectrum and the potential for the case of U = 0, ∆ = 10h 0 .In figure 5, absolute differences between neighboring ∆A i values are shown with vertical dashed lines and we see that these lines overlap nicely with increases in spectral yield.This overlap is well aligned with the qualitative real-space picture of the dynamics introduced in section 3.1.2.As U = 0, the energy transitions associated with nearest-neighbor hopping of the electrons are determined entirely by Ĥ∆ , which means that signatures are expected in the spectrum at these energies.This is so since the hopping of an electron from a site i to j with A i > A j releases an amount of energy corresponding to ∆(A i − A j ).This energy release corresponds to a high harmonic with that given energy.Due to this mechanism we expect, and observe, localized enhancement features in the HHG spectrum centered on the energies corresponding to the differences between the ∆A i for neighboring sites.Here, we simulate 12 sites, i.e. we have significantly fewer  energy gaps in the system than would be the case in an experimental setting.The limited number of sites in the numerical simulation means that the individual energy gaps may be determinable from figure 5, but are expected to be washed out in an experimental setting.

Averages over spectra
To describe the effects of Ĥ∆ on the spectra in a more general sense, we have calculated the spectra for ten different random sets of A i -values, {A i }, each generated according to equation (6).Those spectra are plotted in figure 6.In figure 6(a), U = 0 and ∆ = h 0 , we observe two noteworthy things.Firstly, there are only relatively small differences between the spectra for different {A i } sets, which makes sense as ∆ remains relatively small.Secondly, the general structure of the spectra, particularly the averaged spectra, consists of a plateau that ends at the 20th harmonic order, followed by a dropoff.The bandwidth of the Bloch band plus two times the ∆ value used here corresponds to about 23 harmonic orders.The 20 harmonic orders are, therefore, close to the maximal possible energy transition.As the probability of having picked A i values that facilitate the maximal possible energy transition is rather low, it seems that the energy around 20 harmonic orders is the maximal energy associated with a transition in the lattice for a given random set of {A i }.This observation explains the overall spectral structure.As ∆ increases, see figure 6(b) where U = 0 and ∆ = 5h 0 , the spectral dependence on the A i values increases.This increase is seen from the substantially larger difference between the spectra from different sets of A i -values.This dependence indicates that each spectrum depends significantly on the specific energy gaps for the individual realizations of the A i values and the ability of the pulse to drive electrons past those energy gaps.Notice, in particular, the bump in one of the realizations at around the 30th harmonic order.That feature is similar to those observed in figure 5 and in line with the discussion in section 3.4.Such features will be washed out by averaging over a number of realizations of {A i } values.For larger U, see figures 6(c) and (d) U = 5h 0 and ∆ = h 0 , ∆ = 5h 0 , respectively, we observe less of a dependence on the specific A i values.This lowered dependence illustrates the effect of the ∆ and U terms balancing each other out.In figure 6(e), U = 10h 0 and ∆ = h 0 , we observe the effect of U being much larger than those of ∆, meaning the specific choice of A i values has equivalently little effect.Finally in figure 6(f), U = 10h 0 and ∆ = 5h 0 , the effects of the choice of A i values are again relatively small, which is a result of the comparable ∆ and U values, which means that U can balance out the energy gaps induced by Ĥ∆ .
To explore the sensitivity of the spectra to the method used to simulate the randomness at each site (the A i 's of ( 6)), we have generated spectra similar to those of figure 4 using both a normal distribution and a uniform distribution of A i values, in place of equation (6).We did this with 10 different choices of A i values for each spectrum, and the results, shown and discussed in the appendix, confirm the conclusions drawn above.

Summary, conclusion and outlook
In this work, we used the Hubbard model to study the effects of lattice imperfections on HHG from correlated systems.We simulated lattice imperfections by adding a random modification to the chemical potential at each lattice site, which mirrors the procedure taken in studies of Anderson localization in the uncorrelated limit.Previous studies of HHG from highly correlated systems have, as far as we are aware, operated under the assumption of a perfect lattice.Such studies led to the doublon-holon-based three-step model [21] and similar physical pictures largely based around interaction-induced energy transitions [21,37].As discussed here, lattice imperfections can modify the dynamics of such a system significantly.
We set out to answer the following questions.(i) How are HHG spectra affected by having a lattice with a degree of randomness in the site-dependent potentials in the presence of electron-electron interactions?(ii) How do effects of having such randomness in the lattice and electron-electron interactions relate?And (iii) how do dynamics in the lattice relate to the spectra?We addressed these questions via two qualitative physical pictures of the system: The k-space-based doublon-holon band picture, where the 3-step model of HHG from correlated materials is described using a set of qualitative bands through which doublons and holons propagate.The second real-space-based picture relies on the idea that lattice imperfections lead to ground states with high electron population on a few sites with low chemical potential, while electron-electron interactions tends to drive the system towards a more even population.Based on this picture, we studied the dependence of the Mott gaps on ∆ and the specific set of A i values.We found that the Mott gap is only weakly affected by ∆ when ∆ remains very small but begins to decrease approximately linearly once localization of doublons and holes is achieved.We also found that past ∆ values of ∆ = U/2, the Mott gap becomes more sensitive to the specific choice of A i values and looses meaning as doublons start appearing in the ground state of the half-filled system.
Studying the HHG spectra, we found that the fundamental difference between the effects induced by imperfection and interactions results in an internal struggle between the two effects, resulting in more dynamics and HHG when 2∆ ≈ U.The spectra show significant gain in the medium to high harmonic orders, but less so for the lowest orders, as soon as imperfections are added.This trend continues until the ∆ value becomes very large, at which point the overall spectra begin to lose gain as the electrons become localized to sites with low chemical potential.As electron-electron interaction is added on top of the imperfections the effects of imperfections largely disappears, so long as the interactions remain weak.As U approaches U = 10h 0 , the effects of the imperfections become increasingly noticeable.Severe imperfections result in increased harmonic gain among the lower harmonic, opposite to the effect when the interactions were weak.We linked this to the balancing act described above.
In short, the effects of lattice imperfections on the spectra seem largely similar to those of the correlations, with the interesting add-on of correlations and the imperfections balancing each other resulting in Bloch-like transitions in lattices with very large correlations and imperfections.So if the correlations is low then imperfections result in increased high harmonic gain, whereas if the correlations are large then imperfections result in gain in the lower-order harmonics.
In studies of impurity effects in HHG from systems with significant beyond-mean-field electron-electron interaction, it could be of interest to consider effects of changing the hopping term in additional to the conventional site-dependent changes in chemical potential.Furthermore, one could attempt to change the relative positions of the sites, by introducing site-dependent lattice spacing, or look at the effects of single impurities, by, e.g.modifying the chemical potential at only one location, and potentially neighboring sites.Hence, there are a multitude of possibilities to investigate the mechanisms of HHG from disordered systems with significant beyond mean-field electron-electron interaction.

Appendix. Other distributions for the simulation of imperfections
In this appendix, we explore the effect of sampling the A i values with different statistical distributions than the one presented in equation (6).Specifically, we discuss the equivalents of figure 6, but now generated from a uniform distribution, figure 7, and Gaussian distribution, figure 8, respectively.The uniform distribution of A i -values is limited to the interval [−1, 1], and the Gaussian distribution of the A i values has a standard deviation of 1 and is centered at 0. Note that the centering of the Gaussian has no impact on the dynamics, as it only serves as an energy shift of every state of the system.Starting with the results based on the uniform distribution of figure 7, we observe the same qualitative behavior as observed in figure 6.In figure 7(a), U = 0 and ∆ = 1h 0 , the spectra overlap, which is consistent with the observations presented in the main part of the manuscript, and is expected as ∆ is too small for Ĥ∆ to dominate over ĤHop , which also explains the drop-off at around the bandwidth of the Bloch band.In figure 7(b), U = 0 and ∆ = 5h 0 , there are larger differences between the spectra, which makes sense since Ĥ∆ is relatively more important and therefore the choice of A i values becomes important.In figures 7(c), U = 5h 0 and ∆ = h 0 , (d), U = 5h 0 and ∆ = 5h 0 , and (e), U = 10h 0 and ∆ = h 0 , the spectra largely overlap, because ĤU is sufficiently important to balance out the effects of Ĥ∆ , making the choice of A i values less significant.In figure 7(f), U = 10h 0 and ∆ = 5h 0 , we also observe considerable overlap with the exception of the blue spectrum, which dips under the other spectra from around the 30th harmonic and down.The specific realization of A i belonging to that result has a noticeably more constant set of A i values than the remainder, which makes it distinct from the rest.In conclusion, the results from the uniform distribution of A i values are entirely in line with those presented in the manuscript.
Regarding the Gaussian distribution, it is important to note that the A i values can be larger than 1 or smaller than −1, which they cannot in either of the other sampling methods considered.This difference effectively increases the importance of Ĥ∆ for this distribution of A i values.For comparison, the average variance of the A i values from the uniform distribution is 0.30, whereas from the Gaussian distribution it is 0.88, almost 3 times larger.We observe this increased importance of Ĥ∆ across all the panels of figure 8, as the variance between the spectra from each iteration of A i values are larger in every panel than in the corresponding panels of figures 7 and 6.However, knowing this, one can see that figure 8 is entirely in line with the conclusions presented above.When Ĥ∆ is insignificant compared to ĤHop and ĤU , the choice of A i values has little impact on the spectra, as seen in figures 8(a), U = 0 and ∆ = h 0 , (c), U = 5h 0 and ∆ = h 0 , and (e), U = 10h 0 and ∆ = h 0 .In the opposite case, where ∆ is so large that Ĥ∆ dominates over ĤHop and ĤU , the choice of A i values has a significant impact on the spectra, as observed in figures 8(b), U = 0 and ∆ = 5h 0 , (c), U = 5h 0 and ∆ = 5h 0 , and (f), U = 10h 0 and ∆ = 5h 0 .Note the trend across those panels, (b), (d), and (f), towards increasingly similar spectra, as expected due to the increasing U value.In short, the conclusions presented in the main part of this paper are robust not just to changes in A i values from the distribution presented in the main part, but to different choices of A i values from other distributions too.

Figure 1 .
Figure 1.Histogram of the eigenenergies of the U = 10h0, ∆ = 0 system.The blue arrow indicates the Mott gap, ∆Mott(U, ∆) of equation (10), and the green arrows indicate the width of the Bloch band, ∆ width of equation(11).

Figure 2 .
Figure 2. Illustration of the dominant ground state configurations of systems in the limit of h0 ≪ (U, ∆).The thick black lines represent the modifications to the chemical potential at the various sites, and the thin, short black lines with crosses represent electrons.As illustrated by the double-headed arrow in panel (b), the second electron placed on a given site is displaced upwards by U, in order to represent the energy cost of adding a doublon.

Figure 3 .
Figure 3. Mott gaps calculated for U = 10h0, ∆-values between 0 and 10h0, and 10 different sets of A i -values, see equation (6).Each calculated value of ∆Mott(U, ∆) are highlighted with a star.Each color represents a given set of A i -values.

Figure 5 .
Figure 5. (a) HHG spectrum for the U = 0, ∆ = 10h0 system with the potential shown in (b).The vertical dashed lines in (a) indicate the harmonic orders corresponding to the absolute energy differences between the ∆A i values on each pair of neighboring sites in (b).The energy difference in ∆A i between site 7 and 8 corresponds to harmonic ≃ 73 and is not shown in (a).The thick dashed line at harmonic ≃ 16-17 corresponds three energy differences.

Figure 6 .
Figure 6.Spectra generated from ten different random realizations of A i -values, in accordance with the probability distribution described in connection with equation (6).Such a set of spectra is generated for a variety of U and ∆ values, (a) and (b), U = 0, (c) and (d), U = 5h0, and (e) and (f) U = 10h0.The ∆ values in (a), (c), and (e) are ∆ = 1h0, and in (b), (d), and (f) ∆ = 5h0.The red spectra are the averages of the other spectra.Each set of A i values corresponds to a given color.

Figure 7 .
Figure 7. Spectra generated from ten different random realizations of A i -values, from a uniform probability distribution, limited to the interval [−1, 1].Such a set of spectra is generated for a variety of U and ∆ values, (a) and (b), U = 0, (c) and (d), U = 5h0, and (e) and (f) U = 10h0.The ∆ values in (a), (c), and (e) are ∆ = 1h0, and in (b), (d), and (f) ∆ = 5h0.The red spectra are the averages of the other spectra.Each set of A i values corresponds to a given color.

Figure 8 .
Figure 8. Spectra generated from ten different random realizations of A i -values, from a Gaussian distribution, with a standard deviation of 1.Such a set of spectra is generated for a variety of U and ∆ values, (a) and (b), U = 0, (c) and (d), U = 5h0, and (e) and (f) U = 10h0.The ∆ values in (a), (c), and (e) are ∆ = 1h0, and in (b), (d), and (f) ∆ = 5h0.The red spectra are the averages of the other spectra.Each set of A i values corresponds to a given color.