format: text

text:
# Interaction effects in non-Hermitian models of vortex physics

## I Introduction

An external magnetic field between the lower critical field $`H_{c1}`$ and the upper critical field $`H_{c2}`$ of a type-II superconductor penetrates the material in the form of flexible flux lines (or vortices). In high-temperature cuprate superconductors, which are strongly type-II superconductors, a wide variety of thermodynamic and dynamic phases can be formed in a system of magnetic flux lines by an interplay of strong thermal fluctuations, intervortex interactions, material anisotropy, and disordering and pinning effects generated by various kinds of defects and impurities. In addition to the famous Abrikosov vortex lattice phase, a number of fundamentally new liquid, crystalline and glassy phases have been proposed over the last decade. The study of these phases and the transitions among them is a subject of active theoretical and experimental research in condensed matter physics.

A convenient way of understanding the system of interacting flux lines is provided by the formal mapping between the classical statistical mechanics of ($`d+1`$)-dimensional directed flux lines and the nonrelativistic quantum mechanics of $`d`$-dimensional bosons. In this mapping, the flux lines traversing the sample along the direction of the external magnetic field, $`𝐇=H\widehat{𝐳}`$, correspond to the boson world lines propagating through the imaginary time dimension. The thickness of the sample in the $`z`$ direction, $`L_z`$, is analogous to the inverse temperature $`\beta \mathrm{}`$ ($`=\mathrm{}/k_BT`$), which is the size of the imaginary time dimension, $`L_\tau `$. The effective line tension $`\stackrel{~}{ϵ}_1`$ corresponds to the boson mass $`m`$ and thermal fluctuations proportional to $`k_BT`$ play the role of quantum fluctuations proportional to $`\mathrm{}`$. When the flux lines satisfy periodic boundary conditions in the $`z`$ direction or in the thermodynamic limit $`L_z\mathrm{}`$ (corresponding to the zero temperature limit for the fictitious bosons), the boson-vortex mapping becomes precise. When free boundary conditions are more appropriate to the flux line system, the boson mapping must be corrected by finite size effects. More details of the boson-vortex mapping will be discussed in Sec. II A and summarized in Table I.

Using the mapping between flux lines and bosons, it has been suggested that, in clean high-temperature superconductors, the Abrikosov vortex lattice can melt and be transformed to an entangled flux liquid, which is an analogue of the bosonic superfluid phase. It has also been argued that it is possible to produce new glassy phases by adding various kinds of defects to clean superconductors. In the present work, we focus on the “Bose glass” phase, which can be produced, for instance, by artificially creating parallel columnar defects inside a superconductor using heavy ion irradiation. These defects are assumed to be distributed randomly in the plane perpendicular to the columns. In the spirit of the boson-vortex mapping, the columnar defects provide a static random potential for bosons, if they lie parallel to the external magnetic field. Interacting bosons in a strong static random potential can be localized and form the insulating Bose glass phase. The quantum superfluid to Bose glass transition, occurring when the effective disorder strength is changed across a critical value at zero temperature, has been studied extensively during the last decade. The knowledge obtained from this research can be directly applied to the flux-liquid to Bose-glass transition in high-temperature superconductors with parallel columnar defects.

If the direction of the external magnetic field does not coincide with that of columnar defects, it is convenient to separate the transverse component of the field $`𝐇_{}`$ from the parallel component $`𝐇_{}`$ ($`=H_{}\widehat{𝐳}`$). When $`H_{}H_{}`$, or equivalently when the tilt angle of the external field, $`\theta _H=\mathrm{tan}^1(H_{}/H_{})`$, is sufficiently small, one can show that the $`z`$ axis corresponds to the imaginary time dimension and the transverse component $`H_{}`$ plays the role of an imaginary vector potential for bosons. The corresponding fictitious quantum Hamiltonian is non-Hermitian and in general possesses complex energy eigenvalues.

Recently, there has been great interest in non-Hermitian quantum mechanics of noninteracting bosons in a constant imaginary vector potential, which we call $`𝐠`$, and a random scalar potential. In the boson-vortex analogy, $`𝐠`$ is proportional to $`𝐇_{}`$. Since the Hamiltonian is non-Hermitian, the energy eigenvalues can be either real or complex. Hatano and Nelson have argued that all states with complex eigenvalues are extended, whereas those with real eigenvalues are either localized or extended. It is simpler to describe the basic phenomena in the ($`1+1`$)-dimensional case. A model system (with periodic boundary conditions) appropriate to flux lines in ($`1+1`$) dimensions is sketched in Fig. 1. In the absence of the imaginary vector potential, it is well-known that all eigenstates are localized (in both one and two dimensions) and the corresponding eigenvalues are real. This continues to be the case if $`𝐠`$ is sufficiently weak or the random potential is sufficiently strong. As the imaginary vector potential, or the transverse magnetic field, is increased above a critical value $`g_1`$, there appear a few extended states with complex eigenvalues in the center of the energy band. The spectrum of this group of extended states is shaped like a bubble on a complex energy plane and the two points where the bubble is attached to two real lines of localized states are the mobility edges separating the extended states from the localized ones (see Fig. 2). As $`g`$ increases further, more states change from being localized to extended and the complex bubble expands. Ultimately, above another critical value $`g_2`$ of the imaginary vector potential, all states including the ground state become extended and the entire spectrum is shaped like a single complex bubble. The basic physics is similar in two dimensions, although localized and extended states coexist near the band center and the spectrum is much more complicated.

These results suggest the existence of a novel delocalization transition even in one- and two-dimensional non-Hermitian systems. In a small external transverse magnetic field, the flux lines remain in localized states and are strongly pinned by columnar defects. Since the system displays perfect diamagnetism against the transverse field in this case, this phenomenon has been termed the transverse Meissner effect. However, above a certain critical transverse field, the flux lines are depinned from columnar defects and form a tilted flux liquid. A schematic phase diagram of high-temperature superconductors in the $`T`$-$`H_{}`$ plane is shown in Fig. 3.

In spite of the large number of papers devoted to the non-Hermitian problem, there has been very little research on the effects of interactions on the delocalization phenomena in non-Hermitian systems. One may attempt to take account of the interaction effects by assuming that the nature of the spectrum in the interacting case remains the same as in the noninteracting one and forbidding multiple occupancy of the single-particle localized eigenstates. Let us consider the ($`1+1`$)-dimensional case and assume that the transverse field $`H_{}`$ is in a range such that there is a mobility edge separating low-energy localized states from high-energy extended ones. In the boson-vortex mapping, the parallel component of the field, $`H_{}`$, corresponds to the chemical potential $`\mu `$ of the bosonic system. As $`H_{}`$ is increased at a fixed $`H_{}`$, the states are filled in order of increasing energy to obtain the ground state. At a critical value of $`H_{}`$, all states below the mobility edge are filled and additional flux lines begin to fill the extended states just above the boundary. These extended states describe tilted flux lines and the tilt angle is known to be finite at the mobility edge. The Hartree-Bogoliubov results presented here suggest that this naive picture, which is similar to the “inert layer” picture used to explain the superfluid onset transition in 4He adsorbed in porous media or on random substrates, is wrong and needs to be substantially modified.

Hwa et al. have developed a physical picture of the delocalized phase and the delocalization transition in ($`1+1`$) dimensions in the presence of interactions. When $`H_{}`$ is larger than a critical value, each flux line jumps from one columnar defect to another via kinks. These kinks line up and form chains to reduce the interaction energy cost. The density of the chains goes continuously to zero as $`H_{}`$ decreases to the critical value and the universality class of the delocalization transition is predicted to be the same as that of the standard commensurate-incommensurate transition in ($`1+1`$) dimensions. The physical picture in the ($`2+1`$)-dimensional case is less clear and expected to be more complicated than that in the ($`1+1`$)-dimensional one. More recently, Lehrer and Nelson studied the non-Hermitian Mott-insulator to flux-liquid transition occurring in a system of interacting flux lines with a periodic array of columnar defects, using both the mean-field theory and the renormalization-group method. The generic universality class of the non-Hermitian transition was found to be the same as that of flux lines entering the Meissner phase at the lower critical field $`H_{c1}`$.

The main aim of this work is to develop a quantitative theory for understanding the delocalized phase and the delocalization transition in interacting non-Hermitian systems in the presence of strong disorder. We apply the theory to the disordered flux liquid phase and to the flux-liquid to Bose-glass transition in high-temperature superconductors with columnar defects. Generalizing the model of Ref. 17 to the short-range-interacting case, we introduce a $`d`$-dimensional non-Hermitian Hubbard model for lattice bosons in a random potential:

$``$ $`=`$ $`{\displaystyle \frac{t}{2}}{\displaystyle \underset{𝐱}{}}{\displaystyle \underset{\nu =1}{\overset{d}{}}}\left(e^{𝐠𝐞_\nu /\mathrm{}}a_{𝐱+𝐞_\nu }^{}a_𝐱+e^{𝐠𝐞_\nu /\mathrm{}}a_𝐱^{}a_{𝐱+𝐞_\nu }\right)`$ (2)
$`+{\displaystyle \underset{𝐱}{}}V_𝐱a_𝐱^{}a_𝐱+{\displaystyle \frac{U}{2}}{\displaystyle \underset{𝐱}{}}a_𝐱^{}a_𝐱^{}a_𝐱a_𝐱,`$
where $`a_𝐱^{}`$ and $`a_𝐱`$ are boson creation and annihilation operators at lattice site $`𝐱`$ respectively. The hopping integral $`t`$ sets the scale of the boson “kinetic energy” and $`U`$ is the strength of the repulsive on-site interaction. The random potential $`V_𝐱`$ is assumed to be uncorrelated in space and uniformly distributed between $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }`$. The constant vector $`𝐠`$ is a non-Hermitian external field and the vectors $`\{𝐞_\nu \}`$ are the unit lattice vectors. We assume a hypercubic lattice with the lattice constant $`a`$, therefore, $`|𝐞_\nu |=a`$. It is convenient to introduce a dimensionless non-Hermitian field $`𝐡𝐠a/\mathrm{}`$. The connection between this model and the physics of vortices will be discussed in more detail in Sec. II A. The transverse component of the external magnetic field $`𝐇_{}`$ is related to the non-Hermitian field $`𝐡`$ via

$$𝐡=\frac{\varphi _0𝐇_{}a}{4\pi k_BT},$$
(3)
where $`\varphi _0=\pi \mathrm{}c/e`$ is the magnetic flux quantum. We use boson notation in the discussion below.

The Hermitian version of the boson Hubbard model (2), which is obtained by setting $`𝐠=\mathrm{𝟎}`$, and its continuum counterpart have been studied extensively using a variety of methods, which include numerical simulations, scaling analysis, renormalization group calculations and the Hartree-Bogoliubov approximation.  Inspired by an earlier work by Lee and Gunn, Singh and Rokhsar developed a numerical method of studying the disordered superfluid phase and the quantum superfluid to Bose-glass transition in a system of interacting bosons in a random potential, for the purpose of understanding some experiments done with 4He adsorbed in porous media. Using a combination of the Hartree approximation and the Bogoliubov theory, they calculated the condensate depletion and the superfluid density, as well as the excitation spectrum and the density of states.

In the present work, we generalize Singh and Rokhsar’s method to the non-Hermitian case. The first step is to determine the ground state wave function in a self-consistent manner using the Hartree approximation. We find that the random potential is screened by a (Hartree) repulsive interaction with other particles in the ground state. The screened potential is weaker and smoother than the original random potential. In the next step, we apply the Bogoliubov canonical transformation  to the ground and excited states obtained using the Hartree approximation and derive the Bogoliubov excitation spectrum and the ground state energy. By taking derivatives of the ground state energy $`E_g`$ with respect to the non-Hermitian field, we calculate the quantities such as the imaginary current $`𝐉_I`$ and the superfluid density $`n_s`$.

The imaginary current in the ground state is defined by

$$𝐉_I=\frac{}{𝐠}\left(\frac{E_g}{N}\right),$$
(4)
where $`N`$ is the number of bosons (or flux lines). This quantity corresponds to the transverse component of the magnetic induction, $`𝐁_{}`$, in flux line systems. The vortex tilt angle,

$$\theta _v=\mathrm{tan}^1\left(\frac{B_{}}{B_{}}\right),$$
(5)
is related to $`J_I`$ by

$`\mathrm{tan}\theta _v`$ $`=`$ $`{\displaystyle \frac{}{h}}\left({\displaystyle \frac{E_g}{tN}}\right){\displaystyle \frac{J_I}{J_0}},J_0={\displaystyle \frac{ta}{\mathrm{}}}.`$ (6)
The quantity $`J_0`$ can be expressed as $`\mathrm{}a^1/m`$, since, in the continuum limit, the hopping integral $`t`$ is related to the boson mass $`m`$ (this is the tilt modulus of the vortices) by

$$m=\frac{\mathrm{}^2}{ta^2}.$$
(7)
In a system of volume $`𝒱`$, the superfluid density is defined by

$`n_s`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{ta^2}}\underset{\delta g0}{lim}{\displaystyle \frac{^2}{(\delta g)^2}}\left[{\displaystyle \frac{E_g\left(𝐠+\delta 𝐠\right)}{𝒱}}\right]`$ (8)
$`=`$ $`\underset{\delta h0}{lim}{\displaystyle \frac{^2}{(\delta h)^2}}\left[{\displaystyle \frac{E_g\left(𝐡+\delta 𝐡\right)}{t𝒱}}\right],`$ (9)
where $`\delta 𝐡(a/\mathrm{})\delta 𝐠`$. The direction of the vector $`\delta 𝐠`$ is that of the superflow. In general, the superfluid density depends on the relative angle between the non-Hermitian field $`𝐠`$ and the superflow. The definition (9) agrees precisely with the more conventional definition  of $`n_s`$ based on the dependence of $`E_g`$ on periodic phase variations (of wavevector $`𝐤_0`$) of the wave function:

$`n_s`$ $`=`$ $`{\displaystyle \frac{m}{\mathrm{}^2}}\underset{k_00}{lim}{\displaystyle \frac{^2}{k_0^2}}\left[{\displaystyle \frac{E_g\left(𝐤_0\right)}{𝒱}}\right].`$ (10)
It turns out that the superfluid density is directly proportional to the vortex part of the inverse tilt modulus, $`c_{44}^{v}{}_{}{}^{1}`$, in flux line systems:

$$c_{44}^{v}{}_{}{}^{1}=\frac{n_s}{n^2\stackrel{~}{ϵ}_1},$$
(11)
where $`n`$ is the density of flux lines. In three-dimensional superconductors, $`n_s`$ and $`n`$ represent areal densities.

Both $`𝐉_I`$ and $`n_s`$ are zero in the Bose glass phase and nonzero in the superfluid (or flux liquid) phase. By computing these quantities for various parameter values, we are able to construct a phase diagram associated with the Bose-glass to flux-liquid transition. When combined with the results on the Hartree and Bogoliubov spectra and associated wave functions, these calculations provide a physical picture for the delocalized phase and the delocalization transition substantially different from the noninteracing case, as discussed below.

Let us again consider the one-dimensional problem where there is a mobility edge separating low-energy localized states from high-energy extended ones. As $`H_{}`$ (corresponding to the boson chemical potential) is increased from the Bose glass side, the localized states are filled in order of increasing energy. At the same time, the interaction effect, which is proportional to the density and screens the random potential, becomes stronger. Since the screening effect is most effective for the many-body ground state, this ground state becomes extended one at a critical value of $`H_{}`$, well before the chemical potential hits the mobility edge in the single particle picture, and the system becomes a tilted flux liquid. As $`H_{}`$ is increased further, more low energy states are transformed to extended ones and eventually all low-energy states will become extended. These arguments suggest that, close to the Bose glass phase boundary, the tilted flux liquid is characterized by a low-energy band of Bogoliubov excitations delocalized by hopping, a band of localized excitations, and then a band of more conventional delocalized states. Its spectrum thus contains two mobility edges.

The outline of the paper is as follows. In Sec. II, we develop the basic formalism using the Hartree and Bogoliubov approximations. In Sec. III, we present the numerical results in one dimension on the energy spectra and associated eigenfunctions both in a random potential and in the presence of a single impurity. The nature of the eigenfunctions is studied by computing the participation ratios and the winding numbers. In Sec. IV, we present the numerical results on the vortex tilt angle and the tilt modulus of ($`1+1`$)-dimensional flux line systems both in a random potential and in the presence of a single impurity. We also construct the phase diagram associated with the flux-liquid to Bose-glass (or Bose-insulator) transition. In Sec. V, we conclude the paper with some remarks. In Appendix A, we develop an analytical theory of both continuum and lattice non-Hermitian models of a Bose gas in a weak random potential using the Bogoliubov approximation.

## II Non-Hermitian Hartree-Bogoliubov theory

### A Model

Our main interest is in understanding the physics of flux lines in high-temperature superconductors with columnar defects. Therefore it is useful to clarify the connection between the boson Hubbard model (2) and the model free energy $`F`$ for $`N`$ flux lines in a sample of thickness $`L_z`$ in the $`z`$ direction (perpendicular to the $`\mathrm{CuO}_2`$ planes) in the presence of $`M`$ columnar defects aligned in the $`z`$ direction, but randomly distributed in the $`xy`$ plane:

$`F={\displaystyle \frac{\stackrel{~}{ϵ}_1}{2}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle _0^{L_z}}\left|{\displaystyle \frac{d𝐫_i(z)}{dz}}\right|^2𝑑z`$ (12)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{ij}{}}{\displaystyle _0^{L_z}}𝒰\left(\left|𝐫_i(z)𝐫_j(z)\right|\right)𝑑z`$ (13)
$`+{\displaystyle \underset{i}{}}{\displaystyle _0^{L_z}}V_D\left[𝐫_i(z)\right]𝑑z`$ (14)
$`{\displaystyle \frac{\varphi _0𝐇_{}}{4\pi }}{\displaystyle \underset{i}{}}{\displaystyle _0^{L_z}}{\displaystyle \frac{d𝐫_i(z)}{dz}}𝑑z,`$ (15)
$`V_D\left[𝐫_i(z)\right]={\displaystyle \underset{k=1}{\overset{M}{}}}V_1\left[\left|𝐫_i(z)𝐑_k\right|\right],`$ (16)
where $`𝐫_i`$ and $`𝐑_k`$ are two-dimensional vectors in the $`xy`$ plane. This free energy is a good approximation if the mass anisotropy $`m_z/m_{}`$ in the underlying Ginzburg-Landau model is sufficiently large, as is the case in cuprate superconductors. $`𝒰(r)`$ is the repulsive interaction potential between flux lines, which can be taken to be local in $`z`$, while $`V_1`$ is the (attractive) interaction between a flux line and a columnar defect. The effective line tension

$$\stackrel{~}{ϵ}_1=\frac{m_{}}{m_z}\left(\frac{\varphi _0}{4\pi \lambda }\right)^2\mathrm{ln}\kappa ,$$
(17)
where $`\kappa =\lambda /\xi `$ is the ratio of the London penetration depth $`\lambda `$ and the coherence length $`\xi `$, arises from the tilt energy of the lines.

The canonical partition function for a system of $`N`$ lines is given by the path integral

$$Z=\underset{i=1}{\overset{N}{}}𝒟\left[𝐫_i(z)\right]e^{\beta F\left[\left\{𝐫_i(z)\right\}\right]},$$
(18)
which can be rewritten in terms of the imaginary-time evolution operator $`e^{L_\tau /\mathrm{}}`$ as

$$Z=\psi ^f\left|e^{L_\tau /\mathrm{}}\right|\psi ^i,$$
(19)
where the bra and ket vectors are the initial and final states, respectively. The lattice version of the quantum Hamiltonian appearing in the above equation is precisely the boson Hubbard Hamiltonian (2). The parameters of this model, expressed in the original flux line language, are summarized in Table I.

In the main part of this paper, we will concentrate on the ($`1+1`$)-dimensional case, though the theory developed here can also be applied to ($`2+1`$) dimensions. For convenience, we rewrite the Hamiltonian (2) as

$`\mu N`$ $`=`$ $`{\displaystyle \frac{t}{2}}{\displaystyle \underset{i}{}}\left(e^ha_{i+1}^{}a_i+e^ha_i^{}a_{i+1}\right)`$ (21)
$`+{\displaystyle \underset{i}{}}\left(V_i\mu \right)a_i^{}a_i+{\displaystyle \frac{U}{2}}{\displaystyle \underset{i}{}}a_i^{}a_i^{}a_ia_i,`$
where $`h=ga/\mathrm{}`$. We assume that the total number of lattice sites is $`L`$. The density is then equal to $`nN/L`$.

### B Hartree approximation

In general, right eigenfunctions of a non-Hermitian operator are not simple complex conjugates of the corresponding left eigenfunctions. Let us consider the right Hartree state

$$\psi ^R(r_1,r_2,\mathrm{},r_N)=\varphi _0^R\left(r_1\right)\varphi _0^R\left(r_2\right)\mathrm{}\varphi _0^R\left(r_N\right),$$
(22)
where all bosons are condensed into the same state described by the normalized single-particle right eigenfunction $`\varphi _0^R(r_i)`$, as a variational ground state for (21). The corresponding left Hartree state is

$$\psi ^L(r_1,r_2,\mathrm{},r_N)=\varphi _0^L\left(r_1\right)\varphi _0^L\left(r_2\right)\mathrm{}\varphi _0^L\left(r_N\right).$$
(23)
The many-particle states (22) and (23) are explicitly symmetric under particle exchange.

The expectation value of the Hamiltonian (21) between the left and right Hartree states is

$`\psi ^L|\left(\mu N\right)|\psi ^R=`$ (24)
$`{\displaystyle \frac{tN}{2}}{\displaystyle \underset{i}{}}\left[e^h\varphi _0^L(i+1)\varphi _0^R(i)+e^h\varphi _0^L(i)\varphi _0^R(i+1)\right]`$ (25)
$`+N{\displaystyle \underset{i}{}}\left(V_i\mu \right)\varphi _0^L(i)\varphi _0^R(i)`$ (26)
$`+{\displaystyle \frac{UN^2}{2}}{\displaystyle \underset{i}{}}\left[\varphi _0^L(i)\varphi _0^R(i)\right]^2.`$ (27)
In order to minimize Eq. (27) with respect to the single-particle state $`\varphi _0^L(i)`$, we need to solve the discrete nonlinear Schrödinger equation:

$`{\displaystyle \frac{t}{2}}\left[e^h\varphi _\lambda ^R(i1)+e^h\varphi _\lambda ^R(i+1)\right]+W_i\varphi _\lambda ^R(i)`$ (28)
$`=(\mu +ϵ_\lambda )\varphi _\lambda ^R(i),`$ (29)
where the effective single-particle potential $`W_i`$ is defined by

$$W_i=V_i+UN\varphi _0^L(i)\varphi _0^R(i).$$
(30)
Similarly, we find that the wave function $`\varphi _\lambda ^L(i)`$ satisfies

$`{\displaystyle \frac{t}{2}}\left[e^h\varphi _\lambda ^L(i1)+e^h\varphi _\lambda ^L(i+1)\right]+W_i\varphi _\lambda ^L(i)`$ (31)
$`=(\mu +ϵ_\lambda )\varphi _\lambda ^L(i).`$ (32)
$`\varphi _0^R(i)`$ and $`\varphi _0^L(i)`$ are the ground state wave functions of Eqs. (29) and (32). We assume that the wave functions satisfy the periodic boundary conditions appropriate to Fig. 1

$`\varphi _\lambda ^R(L)`$ $`=`$ $`\varphi _\lambda ^R(0),\varphi _\lambda ^R(L+1)=\varphi _\lambda ^R(1),`$ (33)
$`\varphi _\lambda ^L(L)`$ $`=`$ $`\varphi _\lambda ^L(0),\varphi _\lambda ^L(L+1)=\varphi _\lambda ^L(1).`$ (34)
Due to the simplicity of the Hamiltonian and the boundary conditions, the eigenfunctions $`\varphi _\lambda ^{R,L}(i)`$ and the eigenvalues $`ϵ_\lambda `$ satisfy simple symmetry relationships

$`\varphi _\lambda ^L(i;h)`$ $`=`$ $`\left[\varphi _\lambda ^R(i;h)\right]^{},`$ (35)
$`ϵ_\lambda (h)`$ $`=`$ $`\left[ϵ_\lambda (h)\right]^{}.`$ (36)
The chemical potential $`\mu `$ is chosen so that the (real-valued) eigenvalue $`ϵ_0`$ is equal to zero. This condition ensures that the Bogoliubov quasiparticle spectrum has no gap, in conformity with general theorems. From Eq. (29) or (32), we then find

$`\mu `$ $`=`$ $`{\displaystyle \frac{t}{2}}{\displaystyle \underset{i}{}}\left[e^h\varphi _0^L(i+1)\varphi _0^R(i)+e^h\varphi _0^L(i)\varphi _0^R(i+1)\right]`$ (38)
$`+{\displaystyle \underset{i}{}}V_i\varphi _0^L(i)\varphi _0^R(i)+UN{\displaystyle \underset{i}{}}\left[\varphi _0^L(i)\varphi _0^R(i)\right]^2.`$
The ground state $`\varphi _0^R(i)`$ ($`\varphi _0^L(i)`$) and the excited states $`\{\varphi _{\lambda 0}^R(i)\}`$ ($`\{\varphi _{\lambda 0}^L(i)\}`$) form an orthonormal basis for single-particle states. In Sec. III, we will solve Eqs. (29-32) numerically in a self-consistent manner.

### C Bogoliubov approximation

In the Bogoliubov approximation, one assumes that the ground state (or the condensate) is macroscopically occupied with occupation number $`N_0`$. This assumption is not true in (or sufficiently close to) the Bose glass phase in any dimension. Even though the superfluid to Bose-glass transition is possible in ($`1+1`$) dimensions, phase fluctuations will actually lead to algebraic decay instead of long-range order in the boson order parameter. We illustrate this behavior in Fig. 4, where we plot the occupation number of the Hartree eigenstates, $`N_\lambda =b_\lambda ^Lb_\lambda ^R`$ (see Eq. (41)), measured in the Bogoliubov ground state against the winding number $`N_w`$ of the corresponding Hartree eigenfunctions. We find that the occupation number rises sharply as the ground state with $`N_w=0`$ is approached. Though the ground state is not macroscopically occupied, it still plays a dominant role in long-wavelength physics and the Bogoliubov approximation can be used to calculate properly quantities such as the imaginary current and the superfluid density. Our theory is basically a mean-field approach and is expected to be valid away from the critical region.

From now on, sums over $`\lambda `$ ($`=1,2,\mathrm{},L1`$) will always exclude the condensate. We expand the boson field operator $`a_i`$ in the complete set of right annihilation operators $`b_0^R`$ and $`\left\{b_\lambda ^R\right\}`$:

$$a_i=\varphi _0^R(i)b_0^R+\underset{\lambda }{}\varphi _\lambda ^R(i)b_\lambda ^R.$$
(39)
Similarly, we expand $`a_i^{}`$ in terms of the left annihilation operators $`b_0^L`$ and $`\left\{b_\lambda ^L\right\}`$:

$$a_i^{}=\varphi _0^L(i)b_0^L+\underset{\lambda }{}\varphi _\lambda ^L(i)b_\lambda ^L.$$
(40)
In the Bogoliubov approximation, we replace the operators $`b_0^L`$ and $`b_0^R`$ by $`\sqrt{N_0}`$, where $`N_0`$ is the number of bosons in the condensate. The total number of bosons is the sum of those in the condensate and those not in the condensate:

$$N=N_0+\underset{\lambda }{}b_\lambda ^Lb_\lambda ^R.$$
(41)
Expanding to first order in the depletion of the condensate, we get

$$b_0^L=b_0^R=\sqrt{N_0}=\sqrt{N}\frac{1}{2\sqrt{N}}\underset{\lambda }{}b_\lambda ^Lb_\lambda ^R+\mathrm{}.$$
(42)
Substituting Eqs. (39), (40) and (42) into Eq. (21) and keeping all terms second order in $`b_\lambda ^L`$ and $`b_\lambda ^R`$, we obtain a quadratic Bogoliubov-type Hamiltonian

$`_B\mu N={\displaystyle \frac{UN^2}{2}}{\displaystyle \underset{i}{}}\left[\varphi _0^L(i)\varphi _0^R(i)\right]^2`$ (43)
$`+{\displaystyle \underset{\lambda \lambda ^{}}{}}\left(A_{\lambda \lambda ^{}}b_\lambda ^Lb_\lambda ^{}^R+{\displaystyle \frac{1}{2}}B_{\lambda \lambda ^{}}b_\lambda ^Lb_\lambda ^{}^L+{\displaystyle \frac{1}{2}}C_{\lambda \lambda ^{}}b_\lambda ^Rb_\lambda ^{}^R\right),`$ (44)
where the matrices $`A`$, $`B`$ and $`C`$ are defined by

$`A_{\lambda \lambda ^{}}`$ $`=`$ $`ϵ_\lambda \delta _{\lambda \lambda ^{}}+UN{\displaystyle \underset{i}{}}\varphi _0^L(i)\varphi _0^R(i)\varphi _\lambda ^L(i)\varphi _\lambda ^{}^R(i),`$ (45)
$`B_{\lambda \lambda ^{}}`$ $`=`$ $`UN{\displaystyle \underset{i}{}}\left[\varphi _0^R(i)\right]^2\varphi _\lambda ^L(i)\varphi _\lambda ^{}^L(i),`$ (46)
$`C_{\lambda \lambda ^{}}`$ $`=`$ $`UN{\displaystyle \underset{i}{}}\left[\varphi _0^L(i)\right]^2\varphi _\lambda ^R(i)\varphi _\lambda ^{}^R(i).`$ (47)
The matrices $`B`$ and $`C`$ describe the processes (including effects due to the random potential!) where a pair of particles are scattered out of and into the condensate respectively, while the second term of the matrix A describes single-particle scattering by the condensate.

We diagonalize the Bogoliubov Hamiltonian (44) by introducing a set of quasiparticle creation and annihilation operators $`\gamma _\mu ^{}`$ and $`\gamma _\mu `$ ($`\mu =1,2,\mathrm{},L1`$), defined by the non-unitary canonical transformation

$`\gamma _\mu ^{}`$ $`=`$ $`{\displaystyle \underset{\lambda }{}}\left(u_{\mu \lambda }b_\lambda ^Lv_{\mu \lambda }b_\lambda ^R\right),`$ (48)
$`\gamma _\mu `$ $`=`$ $`{\displaystyle \underset{\lambda }{}}\left(\overline{u}_{\mu \lambda }b_\lambda ^R\overline{v}_{\mu \lambda }b_\lambda ^L\right),`$ (49)
where the coefficients $`\overline{u}_{\mu \lambda }`$ and $`\overline{v}_{\mu \lambda }`$ are defined by

$$\overline{u}_{\mu \lambda }(h)=\left[u_{\mu \lambda }(h)\right]^{},\overline{v}_{\mu \lambda }(h)=\left[v_{\mu \lambda }(h)\right]^{}.$$
(50)
The quasiparticle operators satisfy

$`[_B,\gamma _\mu ^{}]`$ $`=`$ $`\omega _\mu \gamma _\mu ^{},`$ (51)
$`[_B,\gamma _\mu ]`$ $`=`$ $`\omega _\mu \gamma _\mu ,`$ (52)
where $`\omega _\mu `$ is the (complex-valued) quasiparticle excitation energy. The transformation coefficients $`u_{\mu \lambda }`$ and $`v_{\mu \lambda }`$ turn out to obey the eigenvalue equations

$`{\displaystyle \underset{\lambda ^{}}{}}\left(A_{\lambda \lambda ^{}}u_{\mu \lambda ^{}}+B_{\lambda \lambda ^{}}v_{\mu \lambda ^{}}\right)`$ $`=`$ $`\omega _\mu u_{\mu \lambda },`$ (53)
$`{\displaystyle \underset{\lambda ^{}}{}}\left(C_{\lambda \lambda ^{}}u_{\mu \lambda ^{}}+A_{\lambda ^{}\lambda }v_{\mu \lambda ^{}}\right)`$ $`=`$ $`\omega _\mu v_{\mu \lambda }`$ (54)
and the orthonormality conditions

$$\underset{\lambda }{}\left(\overline{u}_{\mu \lambda }u_{\mu ^{}\lambda }\overline{v}_{\mu \lambda }v_{\mu ^{}\lambda }\right)=\delta _{\mu \mu ^{}}.$$
(55)
In Sec. III, we solve these equations numerically to obtain the coefficients $`u_{\mu \lambda }`$ and $`v_{\mu \lambda }`$ and the eigenvalues $`\omega _\mu `$. Finally, using Eq. (38), we find that the diagonalized Bogoliubov Hamiltonian has the form

$`_B`$ $`=`$ $`E_g+{\displaystyle \underset{\mu }{}}\omega _\mu \gamma _\mu ^{}\gamma _\mu ,`$ (56)
$`{\displaystyle \frac{E_g}{N}}`$ $`=`$ $`{\displaystyle \frac{t}{2}}{\displaystyle \underset{i}{}}\left[e^h\varphi _0^L(i+1)\varphi _0^R(i)+e^h\varphi _0^L(i)\varphi _0^R(i+1)\right]`$ (59)
$`+{\displaystyle \underset{i}{}}V_i\varphi _0^L(i)\varphi _0^R(i)+{\displaystyle \frac{Un}{2}}L{\displaystyle \underset{i}{}}\left[\varphi _0^L(i)\varphi _0^R(i)\right]^2`$
$`+{\displaystyle \frac{1}{2n}}{\displaystyle \frac{1}{L}}{\displaystyle \underset{\mu }{}}\left(\omega _\mu A_{\mu \mu }\right),`$
where $`E_g`$ is the ground state energy.

## III Spectra and wave functions

### A Ground state wave function and screening

The first step of our numerical calculation is to solve Eqs. (29-32) with the periodic boundary condition (34). In all of our calculations, we will fix the hopping integral $`t`$ to a constant value. Then there are three natural dimensionless parameters in Eqs. (29-32), that is, the non-Hermitian field $`h`$, the interaction parameter $`Un/t`$ and the disorder parameter $`\mathrm{\Delta }/t`$. It turns out that the effect of disorder is negligible if $`h`$ or $`Un/t`$ is sufficiently large, or of course if $`\mathrm{\Delta }/t`$ is sufficiently small. In these cases, we begin with (normalized) trial condensate wave functions $`\varphi _0^L(i)=\varphi _0^R(i)=1/\sqrt{L}`$, which is uniform in space. Then we solve the nonlinear Schrödinger equations with the screened random potential $`W_i=V_i+Un`$ numerically to find new left and right single-particle eigenstates. We obtain updated trial condensate wavefunctions by mixing the initial guess with the ground state eigenstate of the screened potential, using the so-called Broyden mixing method widely used in self-consistent computations of nonlinear equations. This procedure is repeated until a convergence is achieved.

In the presence of disorder, the ground state wave function is nonuniform in space. We use the solution obtained for a large $`h`$ or $`Un/t`$, or for a small $`\mathrm{\Delta }/t`$ as trial wave functions for the cases with a smaller $`h`$ or $`Un/t`$, or with a larger $`\mathrm{\Delta }/t`$. The convergence is very slow if the system size is large and the effective disorder is strong. For all numerical results presented here, the lattice size is limited to $`L200`$.

In Fig. 5, we plot the local condensate density $`n_0(i)=\varphi _0^L(i)\varphi _0^R(i)`$ in a single realization of the random potential for $`h=0.1`$, $`Un/t=0.5`$, $`\mathrm{\Delta }/t=1`$ and $`L=60`$. We notice that the condensate is indeed highly nonuniform. In fact, the local condensate density is large (small) at the local minima (maxima) of the bare random potential $`V_i`$ as shown in the inset of Fig. 5, so that the screened random potential is smoother than $`V_i`$. In Fig. 6, we compare the screened random potential for the same parameter values as in Fig. 5 with the corresponding bare random potential. The screening is more effective for larger values of $`Un/t`$. We have verified that the screening does not produce long-range correlations in the random potential, as illustrated in Fig. 7.

### B Hartree and Bogoliubov spectra

After calculating the self-consistent eigenvalues and eigenfunctions using the Hartree approximation, we use them to solve the Bogoliubov eigenvalue equations (54) and (55). The eigenvalues $`\omega _\mu `$ obtained from this calculation give the Bogoliubov quasiparticle spectrum.

In Fig. 8(a), we show the Hartree spectra of a 200-site lattice in a single disorder configuration for $`Un/t=0.5`$, $`\mathrm{\Delta }/t=1`$ and $`h=0.8,0.6,0.3`$. For these parameter values, the system is deep inside the superfluid region. When the non-Hermitian field is strong, the spectrum is close to that of a clean system, with no real eigenvalues except at the very bottom and top of the band. All states are extended in this case. As the non-Hermitian field $`h`$ is reduced, some of the eigenvalues near the top of the band become real at some critical value of $`h`$. This band of real-valued eigenvalues corresponding to localized states grows continuously as $`h`$ is further reduced. We note that there is a mobility edge separating low-energy extended states from high-energy localized states.

In Fig. 8(b), we show the Bogoliubov excitation spectra for the same parameter values as in Fig. 8(a). Similarly to the Hartree spectra, there is a low-energy complex bubble corresponding to a band of extended states. For small values of $`h`$, a band of (localized) states with real eigenvalues appears at the top of the bubble. The Bogoliubov spectrum shows a linear behavior near the ground state, similar to results for non-Hermitian “superfluids” without disorder. This linear vanishing of Im $`\omega _\mu `$ with Re $`\omega _\mu `$ is due to the collective phonon excitations characteristic of the superfluid phase. We notice that the effect of disorder is weaker in the Bogoliubov spectra than in the Hartree spectra in the sense that the critical value of $`h`$ where localized states begin to appear at the top of the band is smaller.

Next, we consider the change of the Hartree and Bogoliubov spectra as we cross the Bose glass-superfluid phase boundary by increasing the value of $`h`$, starting from the Bose glass side. When $`h`$ is zero, all eigenvalues are real and, in a sufficiently disordered system, low-lying eigenstates including the ground state are localized. As $`h`$ is increased, a bubble of extended states with complex eigenvalues appears in the center of the band, but the low-energy states with real eigenvalues remain localized. As $`h`$ crosses the phase boundary, the ground state becomes extended and the system is transformed to a superfluid. Within the mean-field theory framework, the condensate is formed, even though it exhibits a slow algebraic decay at $`T=0`$ in ($`1+1`$) dimensions. As $`h`$ is increased further, a few lowest-lying states become extended and form a small bubble in the complex plane. Note that there is a small band of localized excited states sandwiched between two bands of extended states in this parameter regime! The entire spectrum thus contains three mobility edges. In the continuum limit where lattice effects can be neglected, we expect the two low-energy mobility edges survive, while the top one disappears. For larger values of $`h`$, two bubbles grow in size and eventually merge. This “double bubble” scenario is illustrated in Fig. 9 for a single disorder configuration and for $`Un/t=0.05`$, $`\mathrm{\Delta }/t=1`$, $`h=0.3,0.44,0.5`$ and $`L=60`$. Again, the effect of disorder is weaker in the Bogoliubov spectra than in the Hartree spectra, and the characteristic linear behavior in the Bogoliubov spectrum is observed in the small bubble including the ground state. We have checked that this “double bubble” scenario is quite robust and occurs in any realization of the random potential.

### C Participation ratio and winding number

The nature of the Hartree eigenstates can be studied by calculating the participation ratio

$$P_\lambda =\frac{1}{_i\left|\varphi _\lambda ^L(i)\varphi _\lambda ^R(i)\right|^2}$$
(60)
of state $`\lambda `$. This quantity is a measure of the number of sites where the wave function is not negligible. In case of plane waves, $`P_\lambda `$ is equal to the total number of sites and diverges in the thermodynamic limit. For more general extended states, a weaker condition

$$\underset{L\mathrm{}}{lim}\mathrm{ln}P_\lambda \mathrm{ln}L$$
(61)
is applied. For a localized state, $`P_\lambda `$ is proportional to the localization volume and approaches a constant in the thermodynamic limit.

In Fig. 10, we plot the participation ratios for $`Un/t=0.5`$, $`\mathrm{\Delta }/t=1`$, $`h=0.8,0.3`$ and $`L=100,200`$. We find that, for $`h=0.8`$, the participation ratios for all states scale with the size of the system. This can be considered as evidence that all eigenstates in this case are extended. When $`h=0.3`$, the participation ratios for the states having real eigenvalues at the top of the band do not scale with the size of the system, while those for the states with complex eigenvalues and the ground state do. We conclude that the states near the top of the band with real eigenvalues are localized and there exists a mobility edge separating them from low-energy extended states.

A better criterion for determining the nature of Hartree eigenstates in one dimension is provided by the winding number of the Hartree eigenfunction. This winding number is the generalization for non-Hermitian disordered systems of states labelled by wavevectors in translationally invariant problems. When the eigenvalue $`ϵ_\lambda `$ is complex, the eigenfunction $`\varphi _\lambda ^{R,L}(i)`$ rotates an integer times around the origin on the complex plane, as the spatial dimension (which is a circle) is traversed once. This integer rotation or winding number $`N_w`$ is well-defined even in strongly-disordered cases. It has been argued that the states with nonzero winding numbers are always extended.

In Fig. 11(a), we compare the Hartree spectrum and the winding numbers for $`h=0.3`$, $`Un/t=0.5`$, $`\mathrm{\Delta }/t=1`$ and $`L=200`$. In Fig. 11(b), we show the curves for $`h=0.44`$, $`Un/t=0.05`$, $`\mathrm{\Delta }/t=1`$ and $`L=60`$. A bubble of extended states including the ground state is clearly seen. If the eigenvalue is real, the corresponding eigenfunction is also real, and the winding number is undefined. In this case, we arbitrarily set $`N_w=0`$.

### D Non-Hermitian vortex physics with a single impurity

Some features of the Hartree spectrum can be better understood by examining the cases with a single impurity. In Fig. 12(a), we show the spectra in the presence of a single attractive impurity with the strength $`V_d=5t`$. In the noninteracting case, there is a single bound state, which is the ground state. As the interaction parameter $`Un/t`$ increases, this bound state moves toward the band of extended states with higher real parts of the eigenvalues and, at the same time, the highest energy eigenstate becomes a bound state. Above some critical value of $`Un/t`$, the ground state merges into the band of extended states and another bound state appears at the top of the band.

This behavior can be understood from the shape of the screened potential $`W_i`$ shown in Fig. 12(b). As $`Un/t`$ increases, the potential well of the screened potential becomes shallower and the ground state changes from a bound state to an extended one. On the other hand, the inverted screened potential, $`W_i`$, has two potential wells, which can be sufficiently deep to support one or two bound states for intermediate values of $`Un/t`$. This observation explains the appearance of hole bound states at the top of the band. Extending the present argument to the cases with many (but dilute) impurities with random strengths, we can qualitatively understand the Hartree spectra in the presence of a random potential shown in Figs. 8 and 9.

In Fig. 13, we show the spectra and the screened potential in the presence of a single repulsive impurity with the strength $`V_d=5t`$. As $`Un/t`$ increases, the hole bound state with the highest energy moves toward the band of extended states and merges with it at a critical value of $`Un/t`$. This can also be understood from the shape of the screened potential.

## IV Tilt angle, tilt modulus and the Bose-glass to flux-liquid transition

After obtaining the Hartree eigenvalues and eigenfunctions and the Bogoliubov eigenvalues, we use them to compute the ground state energy per particle, $`E_g/N`$, given in Eq. (59). We note that the last term of Eq. (59) has an explicit dependence on the density $`n`$, while the first three terms do not. This suggests that the density is another independent parameter in addition to $`Un/t`$, $`\mathrm{\Delta }/t`$ and $`h`$. For all parameter values used in this work, the last term of Eq. (59) is numerically much smaller than the other terms, unless the density is very small. In the rest of our calculations, we will fix $`n=0.6`$.

At the final stage of our calculations, we take numerical derivatives of the ground state energy with respect to the non-Hermitian field $`h`$ and calculate the vortex tilt angle, $`\theta _v`$, and the vortex component of the inverse tilt modulus, $`c_{44}^{v}{}_{}{}^{1}`$, using Eqs. (6), (9) and (11). In Fig. 14, we plot $`\theta _v`$ and $`c_{44}^{v}{}_{}{}^{1}`$ as functions of $`h`$, which is proportional to the transverse external field $`H_{}`$, for different values of $`Un/t`$. It is clearly seen that there is a transition from the Bose glass phase with $`\theta _v=0`$ and $`c_{44}^v=\mathrm{}`$ to the flux liquid phase with both quantities finite, as $`h`$ is increased through a critical value $`h_c`$. The transition is broader when the interaction is stronger. It has been predicted in Ref. 24 that, close to the phase boundary, $`c_{44}^{v}{}_{}{}^{1}(hh_c)^{1/2}`$. Our calculations are not precise enough to give critical exponents.

In Fig. 15, we show an approximate phase diagram in the $`U`$-$`h`$ plane constructed using the calculations done on a single-size lattice ($`L=30`$) for a single realization of the random potential with $`\mathrm{\Delta }/t=1`$. We chose the first inflection point of the $`c_{44}^{v}{}_{}{}^{1}`$ curve plotted against $`h`$ as the rough phase boundary. As the interaction strength increases, the system becomes effectively less disordered, and therefore the critical value of $`h`$ required to drive the system to the flux liquid phase becomes smaller.

In Fig. 16, we plot the curves for $`\theta _v`$ and $`c_{44}^{v}{}_{}{}^{1}`$ versus $`h`$ obtained for five different realizations of the random potential. Due to the small system size ($`L=30`$) used in the calculations, these quantities show sizable fluctuation effects. The relative fluctuation of the vortex tilt angle is seen to be smaller than that of the inverse tilt modulus.

The Bose-glass to flux-liquid transition can also be driven by changing parameters such as the disorder strength $`\mathrm{\Delta }/t`$ and the density $`n`$. In Fig. 17, we show $`\theta _v`$ and $`c_{44}^{v}{}_{}{}^{1}`$ as functions of $`\mathrm{\Delta }/t`$. A transition is observed at a critical strength of $`\mathrm{\Delta }/t`$. When the non-Hermitian field $`h`$ is bigger, the transition occurs at a higher value of $`\mathrm{\Delta }/t`$. In this figure, we also compare the numerical results with the analytical formulae in the weak disorder limit derived in Appendix A. Since the variance of the random potential, $`V_i^2V_i^2`$, for the particular disorder configuration used in this calculation is approximately 0.419, we relate the disorder parameter in the analytical formulae, $`\mathrm{\Gamma }`$, to the parameter $`\mathrm{\Delta }`$ by $`\mathrm{\Gamma }0.419\mathrm{\Delta }^2`$.

We have also computed the tilt angle and the tilt modulus for the cases with a single attractive impurity. When both the interaction and the transverse external field is very weak, all flux lines can be pinned by one impurity and the system is an insulator. As the interaction strength is increased however, the screening effect causes the lowest-energy bound state to change into an extended one and the system becomes a tilted flux liquid. This is illustrated in Fig. 18 and the corresponding phase diagram is shown in Fig. 19.

## V Conclusion

In the present paper, we have developed a numerical method of studying the strongly-disordered boson Hubbrad model in the presence of a constant imaginary vector potential and a random scalar potential and used it to understand the physics of vortices in high-temperature superconductors with columnar defects. We have found that the interaction causes the screening of a random potential, which is most effective for the ground and low-lying excited states. For sufficiently large values of the effective interaction parameter or the non-Hermitian field, the nature of these states is changed from localized to extended. All physical properties are strongly influenced by these phenomena. In this paper, we presented a detailed study of the vortex tilt angle and the tilt modulus in ($`1+1`$) dimensions.

###### Acknowledgements.

K. Kim is deeply grateful to Jack Lidmar for numerous enlightening discussions. This work has been supported by the NSF through Grant No. DMR97-14725 and through the Harvard MRSEC via Grant No. DMR98-09363 and by the Korea Science and Engineering Foundation through Grant No. 1999-2-11400-005-5.

## A Non-Hermitian Bogoliubov theory of a dilute Bose gas in a weak random potential

In the weak disorder limit, it is possible to obtain analytical results for various physical quantities using the Bogoliubov approximation. In this Appendix, we describe the results for both continuum and lattice models of a dilute Bose gas in a constant imaginary vector potential and a weak random potential. In the section on the continuum model, we define $`𝒱`$ as the volume and $`n`$ as the number of particles per volume. In the section on the lattice model, $`𝒱`$ represents the number of lattice sites and $`n`$ is the number of particles per site. Our conclusions serve as a check on results obtained via coherent state path integrals in Ref. 7. The results for non-Hermitian lattice models with disorder are new.

### 1 Continuum model

We consider a dilute Bose gas of volume $`𝒱`$ in a weak random potential, $`V(𝐱)`$, and a constant non-Hermitian vector potential, $`𝐠`$, in $`d`$ dimensions. The Hamiltonian is given by

$``$ $`=`$ $`{\displaystyle d^dx\psi ^{}\left[\frac{1}{2m}\left(\frac{\mathrm{}}{i}+i𝐠\right)^2+V(𝐱)\right]\psi }`$ (A2)
$`+{\displaystyle \frac{v_0}{2}}{\displaystyle d^dx\psi ^{}\psi ^{}\psi \psi },`$
where $`\psi (𝐱)`$ is the boson field operator and $`v_0`$ is the interaction strength. This model is a continuum version of the lattice model (2). Upon substituting

$`\psi (𝐱)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{𝒱}}}{\displaystyle \underset{𝐤}{}}e^{i𝐤𝐱}a_𝐤,`$ (A3)
$`V(𝐱)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{𝒱}}}{\displaystyle \underset{𝐤}{}}e^{i𝐤𝐱}V_𝐤`$ (A4)
into Eq. (A2), we obtain

$``$ $`=`$ $`{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{\left(\mathrm{}𝐤+i𝐠\right)^2}{2m}}a_𝐤^{}a_𝐤+{\displaystyle \frac{1}{\sqrt{𝒱}}}{\displaystyle \underset{𝐤,𝐪}{}}V_{𝐤𝐪}a_𝐤^{}a_𝐪`$ (A6)
$`+{\displaystyle \frac{v_0}{2𝒱}}{\displaystyle \underset{𝐤,𝐤^{},𝐤^{\prime \prime }}{}}a_𝐤^{}a_𝐤^{}^{}a_{𝐤^{\prime \prime }}a_{𝐤+𝐤^{}𝐤^{\prime \prime }}.`$
We assume that $`V(𝐱)`$ is a Gaussian random potential, characterized by a single parameter $`\gamma `$:

$`\overline{V(𝐱)V(𝐲)}`$ $`=`$ $`\gamma \delta (𝐱𝐲),`$ (A7)
$`\overline{V_𝐤V_𝐪}`$ $`=`$ $`\gamma \delta _{𝐤,𝐪},`$ (A8)
where $`V_{𝐤=\mathrm{𝟎}}=0`$ and the bar indicates averaging over disorder.

Following Bogoliubov, we replace $`a_0`$ and $`a_0^{}`$ by

$$\sqrt{N_0}\sqrt{N}\frac{1}{2\sqrt{N}}\underset{𝐤\mathrm{𝟎}}{}a_𝐤^{}a_𝐤$$
(A9)
and expand the Hamiltonian in terms of $`a_𝐤`$ and $`a_𝐤^{}`$ with $`𝐤\mathrm{𝟎}`$:

$``$ $`=`$ $`𝒱\left({\displaystyle \frac{g^2}{2m}}n+{\displaystyle \frac{1}{2}}v_0n^2\right)`$ (A13)
$`+{\displaystyle \underset{𝐤\mathrm{𝟎}}{}}\left({\displaystyle \frac{\mathrm{}^2k^2}{2m}}+v_0n+i{\displaystyle \frac{\mathrm{}}{m}}𝐤𝐠\right)a_𝐤^{}a_𝐤`$
$`+{\displaystyle \frac{1}{2}}v_0n{\displaystyle \underset{𝐤\mathrm{𝟎}}{}}\left(a_𝐤a_𝐤+a_𝐤^{}a_𝐤^{}\right)`$
$`+\sqrt{n}{\displaystyle \underset{𝐤\mathrm{𝟎}}{}}V_𝐤\left(a_𝐤^{}+a_𝐤\right)+\mathrm{},`$
where $`n=N/𝒱`$. We ignore the cubic and quartic terms represented by $`\mathrm{}`$ and diagonalize the remaining terms by the nonunitary canonical transformation

$`a_𝐤`$ $`=`$ $`{\displaystyle \frac{c_𝐤^R\alpha _kc_𝐤^L}{\sqrt{1\alpha _k^2}}}`$ (A15)
$`\sqrt{n}{\displaystyle \frac{V_𝐤}{\omega _𝐤\omega _𝐤}}\left({\displaystyle \frac{\mathrm{}^2k^2}{2m}}i{\displaystyle \frac{\mathrm{}}{m}}𝐤𝐠\right),`$
$`a_𝐤^{}`$ $`=`$ $`{\displaystyle \frac{c_𝐤^L\alpha _kc_𝐤^R}{\sqrt{1\alpha _k^2}}}`$ (A17)
$`\sqrt{n}{\displaystyle \frac{V_𝐤^{}}{\omega _𝐤\omega _𝐤}}\left({\displaystyle \frac{\mathrm{}^2k^2}{2m}}i{\displaystyle \frac{\mathrm{}}{m}}𝐤𝐠\right),`$
where

$`\alpha _k`$ $`=`$ $`{\displaystyle \frac{1}{v_0n}}\left[{\displaystyle \frac{\mathrm{}^2k^2}{2m}}+v_0n\sqrt{{\displaystyle \frac{\mathrm{}^2k^2}{2m}}\left({\displaystyle \frac{\mathrm{}^2k^2}{2m}}+2v_0n\right)}\right],`$ (A18)
$`\omega _𝐤`$ $`=`$ $`\sqrt{{\displaystyle \frac{\mathrm{}^2k^2}{2m}}\left({\displaystyle \frac{\mathrm{}^2k^2}{2m}}+2v_0n\right)}+i{\displaystyle \frac{\mathrm{}}{m}}𝐤𝐠.`$ (A19)
The left and right annihilation operators $`c_𝐤^L`$ and $`c_𝐤^R`$ satisfy bosonic commutation relations

$$[c_𝐤^R,c_𝐪^L]=\delta _{𝐤,𝐪},[c_𝐤^R,c_𝐪^R]=0,[c_𝐤^L,c_𝐪^L]=0.$$
(A20)
The diagonalized Bogoliubov Hamiltonian has the simple form

$`_B=E_g+{\displaystyle \underset{𝐤\mathrm{𝟎}}{}}\omega _𝐤c_𝐤^Lc_𝐤^R,`$ (A21)
$`{\displaystyle \frac{E_g}{𝒱}}={\displaystyle \frac{g^2}{2m}}n+{\displaystyle \frac{1}{2}}v_0n^2`$ (A22)
$`{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{𝒱}}{\displaystyle \underset{𝐤\mathrm{𝟎}}{}}\left[{\displaystyle \frac{\mathrm{}^2k^2}{2m}}+v_0n\sqrt{{\displaystyle \frac{\mathrm{}^2k^2}{2m}}\left({\displaystyle \frac{\mathrm{}^2k^2}{2m}}+2v_0n\right)}\right]`$ (A23)
$`n{\displaystyle \frac{1}{𝒱}}{\displaystyle \underset{𝐤\mathrm{𝟎}}{}}{\displaystyle \frac{|V_𝐤|^2}{\omega _𝐤\omega _𝐤}}{\displaystyle \frac{\mathrm{}^2k^2}{2m}},`$ (A24)
where $`E_g`$ is the ground state energy and $`\omega _𝐤`$ is the Bogoliubov quasiparticle spectrum. From Eqs. (4), (9) and (A24), we derive the imaginary current,

$`𝐉_I`$ $`=`$ $`{\displaystyle \frac{𝐠}{m}}\{1{\displaystyle }{\displaystyle \frac{d^dk}{(2\pi )^d}}|V_𝐤|^2`$ (A26)
$`\times {\displaystyle \frac{4\mathrm{cos}^2\theta }{\left[\mathrm{}^2k^2/2m+2v_0n+2(g^2/m)\mathrm{cos}^2\theta \right]^2}}\},`$
and the superfluid density,

$`n_s`$ $`=`$ $`n4n{\displaystyle \frac{d^dk}{(2\pi )^d}|V_𝐤|^2\mathrm{cos}^2\varphi }`$ (A29)
$`\times \{{\displaystyle \frac{1}{\left[\mathrm{}^2k^2/2m+2v_0n+2(g^2/m)\mathrm{cos}^2\theta \right]^2}}`$
$`{\displaystyle \frac{8(g^2/m)\mathrm{cos}^2\theta }{\left[\mathrm{}^2k^2/2m+2v_0n+2(g^2/m)\mathrm{cos}^2\theta \right]^3}}\},`$
where $`\theta `$ is the angle between $`𝐤`$ and $`𝐠`$ and $`\varphi `$ is the angle between $`𝐤`$ and $`\delta 𝐠`$. In case of the Gaussian random potential, the integrals in Eqs. (A26) and (A29) can be performed easily. We summarize the results for the disorder-averaged imaginary current $`\overline{𝐉_I}`$, expressed in terms of the dimensionless parameter $`\stackrel{~}{g}=g/\sqrt{mv_0n}`$, below:

$`\overline{𝐉_I}`$ $`=`$ $`{\displaystyle \frac{𝐠}{m}}\left[1\gamma {\displaystyle \frac{1}{2v_0^2n^2}}\sqrt{{\displaystyle \frac{mv_0n}{\mathrm{}^2}}}{\displaystyle \frac{1}{(1+\stackrel{~}{g}^2)^{3/2}}}\right]`$ (A31)
$`\text{if }d=1,`$
$`\overline{𝐉_I}`$ $`=`$ $`{\displaystyle \frac{𝐠}{m}}\left[1\gamma {\displaystyle \frac{m}{2\pi \mathrm{}^2v_0n}}{\displaystyle \frac{2}{\stackrel{~}{g}^2}}\left(1{\displaystyle \frac{1}{\sqrt{1+\stackrel{~}{g}^2}}}\right)\right].`$ (A33)
if $`d=2`$
The disorder-averaged normal fluid density $`\overline{n_n}`$ ($`n\overline{n_s}`$) is listed below:

$`\overline{n_n}`$ $`=`$ $`\gamma {\displaystyle \frac{1}{2v_0^2n}}\sqrt{{\displaystyle \frac{mv_0n}{\mathrm{}^2}}}{\displaystyle \frac{12\stackrel{~}{g}^2}{(1+\stackrel{~}{g}^2)^{5/2}}}\text{if }d=1,`$ (A34)
$`\overline{n_n}`$ $`=`$ $`\gamma {\displaystyle \frac{m}{2\pi \mathrm{}^2v_0}}{\displaystyle \frac{2}{\stackrel{~}{g}^2}}\left[{\displaystyle \frac{1+2\stackrel{~}{g}^2}{(1+\stackrel{~}{g}^2)^{3/2}}}1\right]`$ (A36)
$`\text{if }d=2\text{}\delta 𝐠𝐠,`$
$`\overline{n_n}`$ $`=`$ $`\gamma {\displaystyle \frac{m}{2\pi \mathrm{}^2v_0}}{\displaystyle \frac{2}{\stackrel{~}{g}^2}}\left(1{\displaystyle \frac{1}{\sqrt{1+\stackrel{~}{g}^2}}}\right)`$ (A38)
$`\text{if }d=2\text{}\delta 𝐠𝐠.`$
Eqs. (A34-A38) have been derived previously using a different method.

### 2 Lattice model

We consider the $`d`$-dimensional non-Hermitian Hubbard model for lattice bosons, Eq. (2), in a weak random potential $`V_𝐱`$. Upon substituting

$`a_𝐱`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{𝒱}}}{\displaystyle \underset{𝐤}{}}e^{i𝐤𝐱}a_𝐤,`$ (A39)
$`V_𝐱`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{𝒱}}}{\displaystyle \underset{𝐤}{}}e^{i𝐤𝐱}V_𝐤,`$ (A40)
where $`𝒱`$ is the total number of sites and $`V_{𝐤=\mathrm{𝟎}}=0`$, into Eq. (2), we get

$``$ $`=`$ $`{\displaystyle \underset{𝐤}{}}\left[ϵ_R+iϵ_It{\displaystyle \underset{\nu =1}{\overset{d}{}}}\mathrm{cosh}(g_\nu a/\mathrm{})\right]a_𝐤^{}a_𝐤`$ (A41)
$`+`$ $`{\displaystyle \frac{1}{\sqrt{𝒱}}}{\displaystyle \underset{𝐤,𝐪}{}}V_{𝐤𝐪}a_𝐤^{}a_𝐪`$ (A42)
$`+`$ $`{\displaystyle \frac{U}{2𝒱}}{\displaystyle \underset{𝐤,𝐤^{},𝐤^{\prime \prime }}{}}a_𝐤^{}a_𝐤^{}^{}a_{𝐤^{\prime \prime }}a_{𝐤+𝐤^{}𝐤^{\prime \prime }},`$ (A43)
where

$`ϵ_R`$ $`=`$ $`t{\displaystyle \underset{\nu }{}}\mathrm{cosh}(g_\nu a/\mathrm{})\left[1\mathrm{cos}\left(k_\nu a\right)\right],`$ (A44)
$`ϵ_I`$ $`=`$ $`t{\displaystyle \underset{\nu }{}}\mathrm{sinh}(g_\nu a/\mathrm{})\mathrm{sin}\left(k_\nu a\right).`$ (A45)
Using the Bogoliubov prescription, we obtain an approximate Hamiltonian

$`_B`$ $`=`$ $`𝒱\left[tn{\displaystyle \underset{\nu }{}}\mathrm{cosh}(g_\nu a/\mathrm{})+{\displaystyle \frac{1}{2}}Un^2\right]`$ (A49)
$`+{\displaystyle \underset{𝐤\mathrm{𝟎}}{}}\left(ϵ_R+iϵ_I+Un\right)a_𝐤^{}a_𝐤`$
$`+{\displaystyle \frac{1}{2}}Un{\displaystyle \underset{𝐤\mathrm{𝟎}}{}}\left(a_𝐤a_𝐤+a_𝐤^{}a_𝐤^{}\right)`$
$`+\sqrt{n}{\displaystyle \underset{𝐤\mathrm{𝟎}}{}}V_𝐤\left(a_𝐤^{}+a_𝐤\right),`$
where $`n`$ is the number of particles per site. We diagonalize this Hamiltonian by the canonical transformation

$`a_𝐤`$ $`=`$ $`{\displaystyle \frac{c_𝐤^R\alpha _𝐤c_𝐤^L}{\sqrt{1\alpha _𝐤^2}}}\sqrt{n}{\displaystyle \frac{V_𝐤}{\omega _𝐤\omega _𝐤}}\left(ϵ_Riϵ_I\right),`$ (A50)
$`a_𝐤^{}`$ $`=`$ $`{\displaystyle \frac{c_𝐤^L\alpha _𝐤c_𝐤^R}{\sqrt{1\alpha _𝐤^2}}}\sqrt{n}{\displaystyle \frac{V_𝐤^{}}{\omega _𝐤\omega _𝐤}}\left(ϵ_Riϵ_I\right),`$ (A51)
where

$`\alpha _𝐤`$ $`=`$ $`{\displaystyle \frac{1}{Un}}\left[ϵ_R+Un\sqrt{ϵ_R\left(ϵ_R+2Un\right)}\right],`$ (A52)
$`\omega _𝐤`$ $`=`$ $`\sqrt{ϵ_R\left(ϵ_R+2Un\right)}+iϵ_I.`$ (A53)
The diagonalized Hamiltonian has the form

$`_B`$ $`=`$ $`E_g+{\displaystyle \underset{𝐤\mathrm{𝟎}}{}}\omega _𝐤c_𝐤^Lc_𝐤^R,`$ (A54)
$`{\displaystyle \frac{E_g}{𝒱}}`$ $`=`$ $`tn{\displaystyle \underset{\nu }{}}\mathrm{cosh}(g_\nu a/\mathrm{})+{\displaystyle \frac{1}{2}}Un^2`$ (A57)
$`{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{𝒱}}{\displaystyle \underset{𝐤\mathrm{𝟎}}{}}\left[ϵ_R+Un\sqrt{ϵ_R\left(ϵ_R+2Un\right)}\right]`$
$`n{\displaystyle \frac{1}{𝒱}}{\displaystyle \underset{𝐤\mathrm{𝟎}}{}}{\displaystyle \frac{|V_𝐤|^2}{\omega _𝐤\omega _𝐤}}ϵ_R.`$
From now on, we will restrict our attention to the one-dimensional case. For a Gaussian random potential satisfying

$$\overline{V_𝐱V_𝐲}=\mathrm{\Gamma }\delta _{𝐱,𝐲}$$
(A58)
and in the thermodynamic limit, we obtain

$`{\displaystyle \frac{\overline{𝐉_I}}{J_0}}=\mathrm{sinh}(h)\{1+{\displaystyle \frac{1}{n}}{\displaystyle _0^\pi }{\displaystyle \frac{d\stackrel{~}{k}}{\pi }}\mathrm{sin}^2(\stackrel{~}{k}/2)`$ (A59)
$`{\displaystyle \frac{1}{2n}}{\displaystyle _0^\pi }{\displaystyle \frac{d\stackrel{~}{k}}{\pi }}{\displaystyle \frac{\left[Un+2t\mathrm{cosh}(h)\mathrm{sin}^2(\stackrel{~}{k}/2)\right]\mathrm{sin}(\stackrel{~}{k}/2)}{\sqrt{t\mathrm{cosh}(h)\left[Un+t\mathrm{cosh}(h)\mathrm{sin}^2(\stackrel{~}{k}/2)\right]}}}`$ (A60)
$`\mathrm{\Gamma }{\displaystyle _0^\pi }{\displaystyle \frac{d\stackrel{~}{k}}{\pi }}{\displaystyle \frac{\mathrm{cosh}(2h)+\mathrm{cos}(\stackrel{~}{k})+2}{\left[2Un\mathrm{cosh}(h)+t\mathrm{cosh}(2h)t\mathrm{cos}(\stackrel{~}{k})\right]^2}}\}`$ (A61)
and

$`\overline{n_s}=n\mathrm{cosh}(h)n_{n1}n_{n2},`$ (A63)
$`n_{n1}={\displaystyle _0^\pi }{\displaystyle \frac{d\stackrel{~}{k}}{\pi }}\{\mathrm{sin}^2(\stackrel{~}{k}/2)\mathrm{cosh}(h)`$ (A64)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\left[Un\mathrm{cosh}(h)+2t\mathrm{cosh}(2h)\mathrm{sin}^2(\stackrel{~}{k}/2)\right]\mathrm{sin}(\stackrel{~}{k}/2)}{\left[t\mathrm{cosh}(h)\right]^{1/2}\left[Un+t\mathrm{cosh}(h)\mathrm{sin}^2(\stackrel{~}{k}/2)\right]^{1/2}}}`$ (A65)
$`{\displaystyle \frac{1}{4}}{\displaystyle \frac{t\mathrm{sinh}^2(h)\left[Un+2t\mathrm{cosh}(h)\mathrm{sin}^2(\stackrel{~}{k}/2)\right]^2\mathrm{sin}(\stackrel{~}{k}/2)}{\left[t\mathrm{cosh}(h)\right]^{3/2}\left[Un+t\mathrm{cosh}(h)\mathrm{sin}^2(\stackrel{~}{k}/2)\right]^{3/2}}}\},`$ (A66)
$`{\displaystyle \frac{n_{n2}}{n}}=\mathrm{\Gamma }{\displaystyle _0^\pi }{\displaystyle \frac{d\stackrel{~}{k}}{\pi }}\{{\displaystyle \frac{\mathrm{cosh}(h)\left[3\mathrm{cosh}(2h)+\mathrm{cos}(\stackrel{~}{k})\right]}{\left[2Un\mathrm{cosh}(h)+t\mathrm{cosh}(2h)t\mathrm{cos}(\stackrel{~}{k})\right]^2}}`$ (A67)
$`{\displaystyle \frac{4\mathrm{sinh}^2(h)\left[Un+2t\mathrm{cosh}(h)\right]\left[\mathrm{cosh}(2h)+\mathrm{cos}(\stackrel{~}{k})+2\right]}{\left[2Un\mathrm{cosh}(h)+t\mathrm{cosh}(2h)t\mathrm{cos}(\stackrel{~}{k})\right]^3}}\},`$ (A68)
where $`\stackrel{~}{k}=ka`$. The quantity $`n_{n1}`$ is the contribution to the disorder-averaged normal fluid density due to the breaking of the continuous translational symmetry in a lattice and $`n_{n2}`$ is the contribution due to disorder.

After preforming the integrals, we find

$`{\displaystyle \frac{\overline{𝐉_I}}{J_0}}=\mathrm{sinh}(h)\{1+{\displaystyle \frac{1}{2n}}`$ (A70)
$`{\displaystyle \frac{1}{n\pi }}\mathrm{tan}^1\left[\sqrt{{\displaystyle \frac{\mathrm{cosh}(h)}{\stackrel{~}{u}}}}\right]{\displaystyle \frac{1}{n\pi }}\sqrt{{\displaystyle \frac{\stackrel{~}{u}}{\mathrm{cosh}(h)}}}`$ (A71)
$`{\displaystyle \frac{\mathrm{\Gamma }}{4t^2}}{\displaystyle \frac{2\mathrm{cosh}(h)+2\stackrel{~}{u}+\stackrel{~}{u}\mathrm{s}ech^2(h)}{\left[\stackrel{~}{u}+\mathrm{cosh}(h)\right]^{3/2}\left[\stackrel{~}{u}+\mathrm{sinh}(h)\mathrm{tanh}(h)\right]^{3/2}}}\},`$ (A72)
and

$`n_{n1}={\displaystyle \frac{1}{2}}\mathrm{cosh}(h)`$ (A74)
$`+{\displaystyle \frac{1}{\pi }}{\displaystyle \frac{1}{\mathrm{cosh}^2(h)}}\{\sqrt{\stackrel{~}{u}\mathrm{cosh}(h)}\mathrm{cosh}(2h)`$ (A75)
$`+[\mathrm{cosh}(h)\mathrm{cosh}(2h)\stackrel{~}{u}\mathrm{sinh}^2(h)]\mathrm{tan}^1\left[\sqrt{{\displaystyle \frac{\mathrm{cosh}(h)}{\stackrel{~}{u}}}}\right]\}`$ (A76)
$`{\displaystyle \frac{1}{\pi }}{\displaystyle \frac{\mathrm{tanh}^2(h)}{\stackrel{~}{u}+\mathrm{cosh}(h)}}\{\sqrt{\stackrel{~}{u}\mathrm{cosh}(h)}[{\displaystyle \frac{3}{2}}\stackrel{~}{u}+\mathrm{cosh}(h)]`$ (A77)
$`+[\mathrm{cosh}^2(h)\stackrel{~}{u}^2]\mathrm{tan}^1\left[\sqrt{{\displaystyle \frac{\mathrm{cosh}(h)}{\stackrel{~}{u}}}}\right]\},`$ (A78)
$`{\displaystyle \frac{n_{n2}}{n}}`$ (A79)
$`={\displaystyle \frac{\mathrm{\Gamma }}{t^2}}\{{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{\left[\stackrel{~}{u}+\mathrm{cosh}(h)\right]^{1/2}\left[\stackrel{~}{u}+\mathrm{sinh}(h)\mathrm{tanh}(h)\right]^{1/2}}}`$ (A80)
$`+{\displaystyle \frac{1}{4}}{\displaystyle \frac{\left[2\stackrel{~}{u}\mathrm{cosh}(h)+\mathrm{cosh}(2h)\right]\left[\stackrel{~}{u}\mathrm{cosh}(h)+2\mathrm{cosh}(2h)\right]}{\left[\stackrel{~}{u}+\mathrm{cosh}(h)\right]^{3/2}\left[\stackrel{~}{u}+\mathrm{sinh}(h)\mathrm{tanh}(h)\right]^{3/2}\mathrm{cosh}^2(h)}}`$ (A81)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\left[\stackrel{~}{u}+2\mathrm{cosh}(h)\right]\left[2\stackrel{~}{u}\mathrm{cosh}(h)+\mathrm{cosh}(2h)\right]\mathrm{tanh}^2(h)}{\left[\stackrel{~}{u}+\mathrm{cosh}(h)\right]^{3/2}\left[\stackrel{~}{u}+\mathrm{sinh}(h)\mathrm{tanh}(h)\right]^{3/2}\mathrm{cosh}(h)}}`$ (A82)
$`{\displaystyle \frac{1}{8}}{\displaystyle \frac{\left[\stackrel{~}{u}+2\mathrm{cosh}(h)\right]^2\mathrm{tanh}^2(h)}{\left[\stackrel{~}{u}+\mathrm{cosh}(h)\right]^{5/2}\left[\stackrel{~}{u}+\mathrm{sinh}(h)\mathrm{tanh}(h)\right]^{5/2}\mathrm{cosh}^2(h)}}`$ (A83)
$`\times [8\stackrel{~}{u}^2\mathrm{cosh}^2(h)+8\stackrel{~}{u}\mathrm{cosh}(h)\mathrm{cosh}(2h)`$ (A84)
$`+\mathrm{cosh}(4h)+2]\},`$ (A85)
where $`\stackrel{~}{u}=Un/t`$. In the continuum limit where $`a0`$, $`ta^2\mathrm{}^2/m`$ and $`\mathrm{\Gamma }a\gamma `$, these expressions reduce precisely to Eqs. (A31) and (A34).