Anderson Localization of Composite Excitations in Disordered Optomechanical Arrays

Optomechanical arrays are a promising future platform for studies of transport, many-body dynamics, quantum control and topological effects in systems of coupled photon and phonon modes. We introduce disordered optomechanical arrays, focusing on features of Anderson localization of hybrid photon-phonon excitations. It turns out that these represent a unique disordered system, where basic parameters can be easily controlled by varying the frequency and the amplitude of an external laser field. We show that the two-species setting leads to a non-trivial frequency dependence of the localization length for intermediate laser intensities. This could serve as a convincing evidence of localization in a non-equilibrium dissipative situation.

Introduction: Optomechanics is a rapidly evolving research field at the intersection of condensed matter and quantum optics [1,2]. By exploiting radiation forces, light can be coupled to the mechanical motion of vibration modes. The interplay of light and motion is now being used for a range of applications, from sensitive measurements to quantum communication, while it also turns out to be of significant interest for fundamental studies of quantum physics.
This rapidly developing area has so far mostly exploited the interaction between a single optical mode and a single mechanical mode. Going beyond this, recent theoretical research indicates the substantial promise of so-called optomechanical arrays, where many modes are arranged in a periodic fashion. In such systems, a large variety of new phenomena and applications is predicted to become accessible in the future. These include the quantum many-body dynamics of photons and phonons [3], classical synchronization and nonlinear pattern formation [4][5][6], tunable long-range coupling of phonon modes [7][8][9], photon-phonon polariton bandstructures and transport [10,11], artificial magnetic fields for photons [12], and topological transport of sound and light [13]. A first experimental realization of a larger-scale optomechanical array has recently been presented, involving seven coupled optical microdisks [14]. Even greater potential is expected for implementations based on optomechanical crystals [15][16][17], i.e photonic crystals that can be patterned specifically to generate localized photon and phonon modes.
Given these promising predictions and the rapid experimental progress towards larger arrays, the question of disorder effects now becomes of urgent importance. For example, in the case of optomechanical crystals, experiments indicate fluctuations in the geometry of about 1%, which translate into equally large relative fluctuations of both the mechanical and optical resonance frequencies. This will invariably have a very significant impact on the transport properties. However, gaining a better understanding of disorder effects in the various envisaged applications is only one motivation of the research to be presented here. Of equal, possibly even greater, importance is the opportunity that is offered by optomechanics to create a highly tuneable novel platform for deliberately studying fundamental physical concepts such as Anderson localization [18].
Localization of waves in a random potential is one of the most remarkable and nontrivial interference effects. Initially, it has been studied in electronic disordered systems [19], though this effect applies equally to other types of quantum and even classical waves [20]. By now, localization and related phenomena have been discovered and investigated in photonic systems [21][22][23][24][25][26][27], coupled resonator optical waveguides [28], cold atomic gases [29,30], in propagation of acoustic waves [31] and in Josephson junction chains [32]. Localization can even play a constructive role, namely in random lasing [27,33]. In spite of extensive theoretical efforts, the unambiguous interpretation of experimental manifestations of localization often remains a challenge, even in situations where the ideal version of Anderson localization applies.
Optomechanical arrays enable controlled optical excitation and readout and at the same time promise significant flexibility in their design. However, it is the optical tuneability of the interaction between two different species (photons and phonons) that makes optomechanical systems a unique platform. As we will show in the present Letter, this offers an opportunity to study effects in Anderson localization physics which currently represent a significant challenge even on the theoretical level and will thus open the door towards exploring novel physics that has not been observed so far. Hamiltonian see Refs. [1,34] for details. Heren ν,j ≡ĉ † ν,jĉ ν,j , and c ν,j is the bosonic annihilation operator of either optical, ν = "o", or mechanical, ν = "m", excitations (we set = 1). Due to disorder, the frequencies ω ν,j fluctuate around mean values ω ν,j d = Ω ν . We assume that ω ν,j are independent Gaussian random variables with variances (ω ν,j − Ω ν )(ω ν ,j − Ω ν ) d = σ 2 ν δ j,j δ ν,ν . Eq. (1) is defined in a rotating frame, where the optical frequencies ω o,j are counted off from the laser frequency, ω L [1]. Thus, Ω o indicates the average detuning and can be tuned in situ by varying the laser frequency. The optomechanical couplings g j are proportional to the mean amplitude of the light circulating in the cavity j [1]. Hence, they are also tunable by varying the laser power.
The presence of two-mode squeezing interactions in Eq. (1) can in principle lead to instabilities. We choose Ω o such that these terms are off-resonant and disorder configurations with optical [35] or vibrational instabilities are very rare. We leave their study for a forthcoming paper.
We can describe the full OMA by a Hamiltonian with nearest-neighbor optical, J o , and mechanical, J m , hopping amplitudes: Our model is time-reversal symmetric [36]. Clean polariton bands: In a clean OMA without dissipation (and without squeezing interaction), the photonphonon hybridization leads to a pair of bands with energies whereΩ = (Ω o + Ω m )/2, δΩ = (Ω o − Ω m ) and likewise forJ, δJ. We refer to Ω ± as upper/lower polariton band, respectively. k denotes the wave-vector of polaritonic Bloch states. We focus on the regime where the uncoupled bands overlap, δΩ < 4J. In this case, the polariton bands are separated by a gap if the coupling becomes large enough, g > g min [37], see Fig.1.
Anderson localization of uncoupled excitations: It is well known that even weak disorder leads to a crucial effect in a 1D system: the eigenstates become localized. If g = 0, each subsystem (photon/phonon), is individually described by the 1D Anderson model [18]. The localized states decay exponentially away from their center, o,m are the bare localization lengths for photons and phonons (for g = 0), measured in units of the lattice constant. Using the theory of 1D localization [38], we can approximate the frequency dependence of the localization length: here the dimensionless quantities χ ν ≡ σ ν /J ν and 2 sin[k ν ] are the disorder strength and the bare group velocity, respectively. Eq.(4) is valid for weak (up to moderately strong) disorder [39]. The comparison of Eq.(4) with numerical results is shown below in Fig.3.
In any experiment, localization can be detected if photons and phonons explore the localization length before leaking out, at a rate κ o,m . This holds true if ξ ν < 2J ν | sin(k ν )|/κ ν , which allows us to neglect dissipation in a first approximation [40]. In addition, the sample size L should be larger than the localization length, L max(ξ (0) ν ). For typical L ∼ 100 we need max(ξ (0) ν ) ∼ 10, corresponding to χ ν ∼ 1. Localization in Optomechanical Arrays: At finite photon-phonon coupling, we encounter an Anderson model with two channels. Localization in the symmetric version of this model (with equal parameters of each channel) is well studied and understood [41,42]. However, OMAs do not fall into this universality class since the mechanical band is generically much narrower than the optical one, J m J o . Thus, the hybrid excitations consist of two components with very different velocities. Similar composite quasiparticles are not uncommon, another example is given by cavity polaritons [43,44] including polaritons in a disordered potential [45]. Developing the theory of localization for such non-symmetric systems remains a real challenge, cf. Ref. [46]. The hybrid localized states typically have two localization lengths, ξ 1 < ξ 2 , see the upper panel of Fig.2. For small systems, L < ξ 1 , the excitations do not feel localization and propagate ballistically. Their transmission decays as exp(−L/ξ 1 ) in the range ξ 1 < L < ξ 2 and becomes suppressed as exp(−L/ξ 2 ) at L > ξ 2 . Our numerical analysis shows that the space region where ξ 1 dominates quickly shrinks with increasing g. Therefore, ξ 2 seems more interesting experimentally, and we will focus on this 'large' localization length in the following. We start from a numerical analysis for relatively strong disorder. At the first stage, we neglect disorder-induced fluctuations of g j [47] and use its homogeneous mean value g = const.
The method: The localization length can be obtained, e.g., from the photon-photon transmission, ] is the frequency-resolved retarded Green's function. T oo is defined via the optical power detected on site j at frequency ω L + Ω while a probe laser of the same frequency is impinging on a different site k [48]. For x = |j − k| → ∞, we expect T oo (j, k; Ω) ∝ exp(−2x/ξ 2 ). Thus, the expression for the averaged (inverse) localization length reads We note that the value of ξ 2 (Ω) is the same for other transmission processes (e.g. photon-phonon transmission) [48]. Eq. (5) can be used as a definition even in the presence of dissipation. In the absence of dissipation and instabilities, there is a simpler alternative, namely extracting the localization length directly from the spatial profile of eigenstates [48]. To ensure reliability of results, we have combined both approaches in numerical simulations.
Analysis of numerical results: The upper panel of Fig.2 shows a typical optomechanical eigenstate in the case of small coupling. The excitation frequency has been selected from the tail of the pure mechanical band. Two different slopes, which correspond to two different localization lengths ξ 1,2 , are clearly visible. When g increases and the other parameters of the upper Fig.2 remain unchanged, the region where ξ 1 dominates shrinks [49] and becomes invisible very quickly. In the following, we will concentrate on ξ 2 and will denote it as ξ for the sake of brevity. The lower panel of Fig.2 illustrates the distribution of optomechanical excitations in space and frequency, including the character of excitations (photon vs. phonon).
Here, and in the following, we have displayed numerical results for an illustrative set of parameters: Localization of the optomechanical excitations becomes pronounced at χ o,m ∼ 1. For concreteness, we have chosen equal relative disorder strength, χ o = χ m = 1. In real samples, J o ranges from 1GHz to 10THz (J m : from 100kHz to 1GHz) with the optical disorder being of order 100GHz to 1THz (mechanical: from 10MHz to 100MHz). Thus, our choice of χ o,m falls into the range of experimentally relevant parameters. The optomechanical coupling in our numerics ranges from weak, g = 10 −3 Ω m , to strong, g = 0.05Ω m . To suppress finite size effects, we employed large systems, L = 10 3 ξ, during exact diagonalization. The Green's functions method has allowed us to explore even much larger sizes.
In Fig.3, we display the frequency-dependence of the localization length of hybrid optomechanical excitations in a disordered array, one of the central numerical results of this article. For comparison, we also show the situation for the uncoupled systems, including the (scaled) analytical expression for ξ  Fig.3a). Once the subsystems are coupled, significant changes of ξ(Ω) occur in the vicinity of the unperturbed (narrow) mechanical band where the optomechanical hybridization is most efficient. Firstly we note that, if loc , the coupling between the optical and the mechanical systems is perturbatively weak even in the middle of the mechanical band [region I in Fig.3(f)]. On the other hand, when the optomechanical coupling  Analytical methods which would allow one to explore localization in strongly disordered systems are not available in general. Nevertheless, it turns out that our optomechanical array corresponds to a certain two-channel system, which was studied analytically in Ref. [46] for the limit of weak disorder and large coupling. Remarkably, the shape of our numerically extracted ξ(Ω) at large g agrees with the predictions of Ref. [46], even though we are here dealing with strong disorder, χ ν ∼ 1 [42]. The theory of Ref. [46] is valid if g is large compared with the (bare) mean level spacing in the localization volume, ∆ (ν) loc , which holds true for the parameters of our numerical study at g ≥ 0.05Ω m [51]. If g > g min , (i.e., if the clean polariton bands are separated by the gap of the width Ω + (k = 0) − Ω − (k = π)) we can use the following (leading in χ ν ) expression for the localization length [46]: Here χ = χ o /C = χ m /C and k ± (Ω) denotes the inverted dispersion relation Ω ± (k). The quantity V ± ≡ 2 sin [k ± (Ω)] is called "rapidity". It coincides with the group velocity of the excitations for g = 0, and according to Eq.(6) it governs the frequency-dependence of ξ(Ω) in the coupled case. The factor C reflects renormalization of the disorder strength caused by the optomechanical coupling. Calculation of C is beyond the scope of Ref. [46] and we have found its approximate value C 1.16 by fitting the analytically calculated maximal value of ξ(Ω > Ω m ) to the numerical one. Fig.3e shows the comparison of the analytical and numerical results. They differ noticeably only close to edges of the clean band where the analytical theory looses its validity because ξ → 1. In addition, the gap is smeared by the relatively strong disorder.
We have discovered that, at Ω Ω m , the crossover between small and large values of g is highly non-trivial (and it is outside the scope of the analytical theory): when the optomechanical coupling increases from g ∼ ∆ This non-trivial dependence of the localization length on the coupling constant, i.e., on the tuneable intensity of the external laser, could help to distinguish localization and trivial dissipation effects in real experiments.
We have checked that the shape of ξ(Ω) is robust with respect to dissipation effects as long as the mean level spacing in the localization volume of the hybrid excitations is larger than the optical and mechanical decay rates, κ ν [53]. Propagation of the excitations is suppressed due to their finite life time which is reflected by the frequency-independent decrease of ξ. Typical profiles ξ(Ω) are shown in Fig.4 where the optical dissipation rate increases until κ o = 0.1Ω m . These profiles are also robust with respect to the spatial inhomogeneity of g j which results from randomness of the cell frequencies [47,54]. important parameters can be easily fine-tuned. Thus, OMAs provide a unique opportunity to study Anderson localization of composite particles in real experiments. Moreover, they should allow to reliably distinguish localization from trivial dissipation effects. Future studies may address the additional novel physics that will arise when two-mode squeezing processes become relevant. At strong driving, this could involve the interplay between instabilities and localization, with interesting connections to random lasing, extending the new research domain of disordered optomechanical arrays into the nonlinear regime. We acknowledge support from the EU Research Council through the grant EU-ERC OPTOMECH 278320. TFR acknowledges support from FAPESP. We are grateful to Vladimir Kravtsov and Igor Yurkevich for useful discussions.
to suppression of the region where ξ1 dominates because the eigenstate is normalized.
[50] The accuracy of Eq.(4) is insufficient to reproduce the maximal value of ξ (0) o (Ωo) at the relatively strong optical disorder, χo = 1. Therefore, we have scaled the analytical answer by a constant factor to adjust the heights.
[51] The mean level spacing in the localization volume is obtained from the relation ∆ loc /∆ = N/ξ . Using the estimate for the mean level spacing in a band of extremely localized states ∆ ∼ δ/N [55, 56], we find ∆ denotes bare (at g = 0) localization lengths. For the chosen parameters, this yields the condition g > 10 −1 max{Jν } ∼ 10 −2 Ωm. Thus, the analytical theory from Ref. [46] can be used to understand the case g = 0.05Ωm while smaller values of g are beyond the validity range of the analytic expressions.
[53] Equivalently, we can require that the localization length is much smaller than the escape length.

Standard optomechanical Hamiltonian
The linearized Hamiltonian, Eq.(2) of the main text, is derived starting from the Hamiltonians of a phononic and a photonic array. In the tight-binding approximation, both Hamiltonians take the same form, Here, the index ν = "o" and ν = "m" refers to the optical and mechanical degrees of freedom. The operatorsĉ ν,j denote the annihilation operator for the site j; ω ν,j and J ν denote random on-site frequencies and constant overlap integrals, respectively. The optical and mechanical modes co-localized on the same site are coupled by the radiation pressure force. The resulting interaction readŝ where g 0 is the eigenfrequency shift of a localized optical mode by a single phonon on the same site. The presence of a laser drive of frequency ω L is described by the additional Hamiltonian term The dynamics of the OMA with the HamiltonianĤ =Ĥ o +Ĥ m +Ĥ om +Ĥ laser is most conveniently described in the rotating frame defined by the unitary transformation We decompose the operatorsĉ ν,j as sums of their mean filed values and new displaced operators, δĉ ν,j , incorporating the fluctuations,ĉ ν,j = ĉ ν,j + δĉ ν,j . After inserting this decomposition intoĤ, all linear terms in δĉ ν,j cancel out in the Hamiltonian. We, thus, reproduce Eq. (7) Here g j ≡ g 0 ĉ o,j are the couplings of the linearized interaction. In the main text, we have investigated a parameter regime where the laser is red-detuned compared to all optical resonances. We have also focused on the OMAs where the broadening of the resonances (set by the typical decay rate κ o ) is smaller than the minimal detuning. In this case, all g j are real valued, consequently, the time-reversal symmetry is preserved by the OM interaction. Summing all contributions, we obtain the HamiltonianĤ of the main text. There, for brevity, we useĉ ν,j and ω o,j for the fluctuation operators and the detunings, respectively. We follow this convention also below. In the main text, we have also assumed that all linearized couplings g j are approximately equal, g j ≈ g. This approximation holds when the mean value of the onsite detuning is much larger than its typical fluctuations, the optical hopping rate, and the typical optical decay rate. Below, we go beyond this approach investigating fluctuating coupling constants g j .

Input/Output formalism
Let us express the elastic part of the photon-photon transmission T oo (k, j; Ω) in terms of the retarded (photonphoton) Green's function.
The response of the OMA to an additional probe field is described by the standard Langevin equations [S1]: whereĉ in ν,j is the input (probe) field. The corresponding output fieldĉ out ν,j is given by the input/output relationŝ For a probe laser of frequency ω L + Ω (corresponding to the frequency Ω in the rotating frame) applied at site k we haveĉ where α p is the amplitude of the probe laser. The transmission is defined with the help of the ratiô The Hamiltonian has been linearized, hence, the response ofĉ o,j to the input field in Eq. (12) is linear. Moreover, for the purpose of calculatingĉ o,j , we can replace the operatorsĉ in ν,j in the source terms of the Langevin equation (12) with their mean valuesĉ in ν,j . We can even formally replace all the input terms with the coherent interaction Thus,ĉ o,j is given by the Kubo formula where H I plays the role of the perturbation. We note that the optomechanical coupling does not conserve the number of excitations and, therefore,ĉ o,j has both an elastic (frequency Ω in the rotating frame or Ω + ω L in the laboratory frame) and an inelastic (frequency −Ω in the rotating frame or −Ω + ω L in the laboratory frame) components. The elastic part of transmission is obtained after time averaging: Using Eqs. (15)(16) and the Kubo formula, we find For j = k, we recover the formula for T oo (j, k; Ω) given in the main text. If the eigenstates of the OMA are localized then T oo (j, k; Ω) decays exponentially on large distances and the inverse localization length can defined as follows: We note that Eq. (19) contains only matrix elements of the Green's function relating operators from sites j and k. These elements can be calculated iteratively with the help of the Dyson's equation which is similar to that suggested in Ref. [S2] for the transfer matrix (details of the algorithm can be found in Ref. [S3]). Generically, one can introduce 4 transmissions in the elastic channel, T νν (j, k; Ω), and 4 transmissions in the inelastic one,T νν (j, k; Ω) (e.g., the photon-photon transmission, inelastic photon-photon transmission). Similar to Eq.(18), these transmissions are described by entries of the matrix Green's function constructed from four-component "Opto-mechanical×Nambu"-spinorsĈ j : parameters of the OMA given in the main text and compared 16 largest localization lengths governed by inserting each component ofĜ R into Eq. (19). These localization lengths coincide up to small numerical errors ≤ 1%, see blue dots in Fig.5. This is the related to the symmetries, namely, particle-antiparticle symmetry, time-reversal symmetry, and space-inversion symmetry. The latter appears effectively in the long disordered OMAs due to the self-averaging. The equivalence of the different transmissions on large distances allows one to find the largest localization lengths of the OMA from any convenient linear combination of |Ĝ R ab (j, j ; Ω)| ensuring a good convergence of the numerical algorithm. In particular, we can use "the generalized transmission" of the opto-mechanical excitations and, after disorder averaging, arrive at: Eq. (23) has been used in the numerical code with the disorder averaging being substituted by the self-averaging of ξ in very long systems, see red dots in Fig.5.

Bogoluibov eigenstates
In the absence of dissipation, there is a simple method which allows one to find the localization length directly from the average number of the excitation. This approach is realized after diagonalizing the Hamiltonian (or, equally, the Heisenberg equations of motion) with the help of the Bogoliubov transformation. Let us define eigenmode operatorŝ d s . Generically,d s can be written as follows: Here u j,s are the Bogoliubov coefficients. The transformation matrix that diagonalize the Heisenberg equations reads as: where U (ν) and V (ν) are N × 2N matrices whose entries are the coefficients u (ν) j,s and v (ν) j,s from Eq. (24). In the absence of dissipation and instabilities, these coefficients satisfy the following relation: The minus sign in front of summands v is caused by the bosonic commutation relations of the operatorŝ c ν,j . As a consequence, T is pseudounitary with the inverse matrix The time evolution of the operatorsĉ ν,j can be obtained using Eq. (27), and it is given bŷ Here ε s is the frequency of the hybrid (opto-mechanical) eigenmode s. Using Eqs. (28), one can derive and find the total average number of excitations at a given site j after the eigenmode s is excited: We have subtracted the (background) fluctuations ofn j in the ground state since the number of excitation in the OMA always fluctuates due to the optomechanical coupling, see the second term in the RHS of Eq.(1). The eigenmodes of the OMA can be found via the numerical diagonalization of the Heisenberg equations. Now we substitute n j (s) for T oo (j, k, Ω) in Eq.(5) and associate the frequency Ω with ε s and the origin k with the coordinate where the eigenmode s has maximal amplitude. This yields the second expression for ξ 2 . Thus, the localization length can be estimated from a log-linear fit of n j (s), see the discussion of Eq.(5) in the main text. Let us analyze the behavior of ξ 2 (g) for the case g ≤ J m √ J m J o where the influence of the optomechanical coupling on the band structure is negligible. The numerical analysis shows that ξ 2 (g) is sub-linear, see the left panel of Fig.6, which indicates the presence of non-perturbative contributions. The full theory for this is missing and we give only phenomenological arguments which are similar to those of the Mott theory [S4] and allow one to explain the sub-linear growth of ξ 2 when g increases up to J m . For simplicity, we concentrate on the transmission T om in the regime ∆ (m) Other parameters correspond to Fig.3a in the main text. Finite transmission T om requires hybridization of bare optical and mechanical states. The main idea of the phenomenological approach is to find a pair of the optical-and the mechanical-states which, being strongly hybridized, provides the largest possible increase of ξ 2 . In other words, we have to estimate the maximal distance between bare localization centers which does not violate the necessary condition for the strong hybridization.
Consider an optical state with the frequency inside the unperturbed mechanical band, Ω m − J m < o < Ω m + J m , see the blue wave-function in the right panel of Fig.6. The space coordinates will be counted from the localization center of this optical state. Such a state can be strongly hybridized with the mechanical states if inequality holds true. Here j is the number of the mechanical state with the localization center at x j > 0 and with the frequency m,j ; M j |O is the overlap between the localized optical and the localized mechanical states Loc. length ξ X 1 2 We recall that ξ for Ω o , cf. Fig.3a. Firstly we note, that, unlike the Mott theory, frequencies o and m,j are not correlated at g = 0. Therefore, | o − m,j | can be arbitrary small even if the localization centers of the bare states are close to each other, x j ξ (0) o,m ⇒ M j |O ∼ 1, see the orange state No.1 in the right panel of Fig.6. On the other hand, it is clear that the 1st mechanical state is unable to support an essential increase of the transmission T om beyond the bare localization length.
T om can become more long-ranged if ξ  Fig.6. However, there is always an optimal state for which x j is relatively large and the smallness of M j |O in Eq.(31) is compensated by the smallness of the frequency separation: cf. the orange state No.3 in the right panel of Fig.6. Now we can speculate that, if ∆ (m) loc ∼ J m /ξ (0) m g J m , T om on large distances and, correspondingly, ξ 2 are governed by the "double-hump" optomechanical state originating mainly from hybridization of the optical state with the optimal mechanical one, see an example in the right panel of Fig.6. Therefore, the largest localization length can be estimated as (C 1,2 are constants of order O(1) which cannot be determined in the frame of the phenomenological approach). The optomechnical state may be "multiple-hump" if x 3 covers several localization volumes of the bare mechanical states. If g ∆ (m) loc we expect a crossover to the purely perturbative regime. Since ∆ (o) loc J m ∼ σ m , the minimal space distance between two optical states belonging to the frequency range of the mechanical band is large. We estimate it as ξ m . This is not the case for the parameters of the main text, in particular, because the disorder is strong.
The fitting in the left panel of Fig.6 can be done equally by using either Eq.(34) or a power-law dependence with some non-universal exponent α < 1. More rigorous theory of the weak coupling regime can be developed by exploiting basic ideas of the virial expansion, see Refs.

Localization of hybrid excitations in the case of fluctuating coupling constant
In the main text, we have concentrated on the case where the optomechanical coupling g is one and the same for all cells. In reality, the coupling fluctuates: g j depends on the mean occupation number of the photons on the site j, g j ∝ c o,j , while the cells with smaller optical frequencies host more photons, see Fig.7. This locally enhances g j on these sites. One can speculate that the coupling constant acquires an effective frequency dependence; g becomes larger for smaller frequencies and it slightly decreases with increasing the frequency. We note that α(ω 0,j ) depicted in Fig.7 is defined as α(ω 0,j ) = c o,j | g0=0 . We do not distinguish α(ω 0,j ) and c o,j since their difference is small, (α(ω 0,j ) − c o,j ) ∼ O(g 2 0 ).  We have recalculated curves from Fig.4 for fluctuating g j . Results are shown in Fig.8. The disorder averaged coupling has been adjusted to the same values (a) g d = g = 0.001Ω m ; and (b) g d = g = 0.05Ω m . Solid/dashed lines show the frequency dependence of the localization length with/without fluctuations of the coupling. The profiles at smaller g are almost intact by its fluctuations. When the mean coupling is larger, the left (right) maximum of ξ 2 (Ω) is suppressed (enhanced) by these fluctuations. This can be explained if, based on Fig.7, we assume that g becomes effectively larger (smaller) at Ω < Ω (Ω > Ω) and notice that the peaks decrease when g increases, cf. Fig.3. Thus, the fluctuations of g j are able to modify the shape of ξ 2 (Ω) at relatively large values of the mean coupling constant g though all qualitatively important features are expected to be robust.