The effects of disorder in dimerized quantum magnets in mean field approximations

We study theoretically the effects of disorder on Bose-Einstein condensates (BEC) of bosonic triplon quasiparticles in doped dimerized quantum magnets. The condensation occurs in a strong enough magnetic field Hc, where the concentration of bosons in the random potential is sufficient to form the condensate. The effect of doping is partly modeled by delta - correlated disorder potential, which (i) leads to the uniform renormalization of the system parameters and (ii) produces disorder in the system with renormalized parameters. These approaches can explain qualitatively the available magnetization data in the Tl_(1-x)K_(x)CuCl_3 compound taken as an example. In addition to the magnetization, we found that the speed of the Bogoliubov mode has a peak as a function of doping parameter, x. No evidence of the pure Bose glass phase has been obtained in the BEC regime.


Introduction
The effects of disorder on the properties of Bose-Einstein condensates present interesting problems both for theoretical and experimental physics [1,2,3,4,5,6]. Disorder is important in various systems of real particles such as superfluid 4 He, cold atoms in optical lattices, and quasiparticles such as polaritons [7] and excitons [8]. These systems are well-suited for experimental studies, however the theory of disordered ensembles of interacting bosons is complex and there are essentially no exact solutions even in one dimension [9]. To approach this problem, Yukalov and Graham (YG) developed a self-consisted stochastic mean field approximation (MFA) [10] for Bose systems with arbitrarily strong interparticle repulsion and arbitrary strength of disorder potential. It was shown that, in general, the Bose system consists of following coexisting components: the condensate fraction, ρ 0 , the normal fraction ρ N , the glassy fraction ρ G , and, in addition, can be characterized by the superfluid density ρ s . In the limit of asymptotically weak interactions and disorder the known results, obtained in pioneering work by Huang and Meng [11] (HM) are reproduced by the YG theory. An interesting question here concerns the problem about the existence of a pure Bose glass (BG) phase, i.e. the phase where the condensate fraction is nonzero, while the superfluid fraction is not yet present. Note that Ref. [6] introduced an alternative definition of the gapless BG phase, having localized short-lived excitations and vanishing superfluid density with a continuous transition to the normal phase at finite temperature. Even without disorder, the condensate is depleted by particle-particle interactions and temperature. The inclusion of random fields depletes the condensate further and, possibly, creates the glassy fraction.
As it was understood recently, a new class of BEC can be provided by spinrelated quasiparticles in magnetic solids such as intensively pumped magnons [12] or triplons in the dimerized quantum magnets in the equilibrium [13]. In the magnets, the effect of disorder, which can be produced by admixing other chemical elements, can be rather strong to be seen in the physical properties such as the temperature-dependent magnetization. The so far most investigated compound showing BEC of triplons is TlCuCl 3 . To study the effect of disorder, solid solutions of quantum antiferromagnets TlCuCl 3 and KCuCl 3 , i.e. Tl 1−x K x CuCl 3 have been experimentally investigated recently [14,15,16] at low temperatures T . The zero-field ground states of TlCuCl 3 and KCuCl 3 are spin singlets with excitation gaps ∆ st = 7.1 K and ∆ st = 31.2 K, respectively and the magnetic excitations are spin triplets. Triplons arise in magnetic fields H > H c , where H c is defined by condition of closing the gap by the Zeeman splitting, that is ∆ st = gµ B H c , where g is the electron Landé factor and µ B is the Bohr magneton. In the mixture Tl 1−x K x CuCl 3 the induced magnetization M exhibits a cusplike minimum at a critical temperature T c (H) for fixed magnetic field H ≥ H c similarly to the parent compound, and can be successfully explained in terms of triplon BEC [17,18,19].
For a theoretical description it is natural to assume that for weak doping x ≪ 1 in the mixed system Tl 1−x K x CuCl 3 a small admixture of potassium forms a disorder potential. Consequently, the recently developed theories of "dirty bosons" [1,2,10,11] can be applied to study the BEC of triplons in Tl 1−x K x CuCl 3 . Here the following natural questions arise. For example, what is the correspondence between admixing parameter x and the properties of the disorder potential ? What are the experimental consequences of the disorder ? Yamada et al. [15] analyzed the electron spin resonance spectrum in Tl 1−x K x CuCl 3 and concluded that there is a Bose glass -BEC transition near a critical magnetic field. Although this interpretation might need a further analysis (see discussion in Refs. [20]) it would be interesting to study the influence of the glassy phase, or more exactly, of the glassy fraction ρ G on the magnetization. Note that even the existence of a pure Bose glass phase still is a matter of debate even in theoretical approaches. For example, it may be predicted by the approach used by Huang and Meng [11] if one extends their formulas from weak disorder to a strong one. On the other hand, no pure Bose glass was found in Monte -Carlo simulations [21] for atomic gases, but predicted for triplons at T = 0 by Nohadani et al. [22].
Here we develop a theory of the disorder effects on the BEC of triplons taking Tl 1−x K x CuCl 3 as a prototype for studies of specific properties. For example, in atomic gases considered in Refs. [10,11], the chemical potential µ is determined self-consistently with fixed number of atoms, while in the triplon gas the chemical potential is a given external parameter controlled by the applied magnetic field and the number of triplons is conserved in the thermodynamic limit. To clarify the terms, we underline that the number of magnons may vary but that of triplons may be tuned and kept fixed, which makes possible the BEC of the latter.
The paper is organized as follows. In Sections II and III we outline the YG and HM approaches valid only for T ≤ T c and extend it for the triplon system. The shift of T c due to disorder and the normal phase properties will be discussed in Section IV. Our numerical results will be presented in Section V. Conclusions will summarize the results of this work.

Yukalov-Graham approximation for disordered triplons
In the following we reformulate the Yukalov-Graham approximation to the triplon system with arbitrary disorder. The Hamiltonian operator of triplons with contact interaction and implemented disorder potential V (r) is given by where ψ(r) is the bosonic field operator, U is the interparticle interaction strength, and K is the kinetic energy operator which defines the bare triplon dispersion ε k . Since the triplon BEC occurs in solids, we perform integration over the unit cell of the crystal with the corresponding momenta defined in the first Brillouin zone. Below we assume that the bare spectrum remains coherent in the presence of disorder and consider it as a simple isotropic one: ε k = k 2 /2m, where m is the triplon effective mass. The distribution of random fields is assumed to be zero -centered, V (r) = 0, and the correlation function Here and below we adopt the units k B ≡ 1, ≡ 1, and V ≡ 1 if not stated otherwise for the unit cell volume.
To describe Bose condensed system where the global gauge symmetry is broken, one employs the Bogoliubov shift: where the condensate density ρ 0 (r) is constant for the homogeneous system, ρ 0 (r) ≡ ρ 0 . Since by the definition the average of ψ † (r)ψ(r) is the total number of particles: with the density of triplons per unit cell ρ = N/V , from the normalization condition one immediately obtains Therefore the field operator ψ 1 (r) determines the density of uncondensed particles. The YG approximation is formulated in representative ensemble formalism, which includes two Lagrange multipliers, µ 0 and µ 1 , defined as: where Ω is the grand thermodynamic potential. It was shown that disorder would not change the explicit expressions for chemical potentials, obtained earlier [23] in Hartree-Fock-Bogoliubov (HFB) approximation without disorder, where σ = 1 V V d 3 r ψ 1 (r)ψ 1 (r) is the anomalous density. The total system chemical potential mu, related to the total number of particles as N = −∂Ω/∂µ, is determined by Clearly, when the gauge symmetry is not broken, i.e. ρ 0 = 0, σ = 0, ρ 1 = ρ, both µ 0 and µ 1 coincide giving µ = µ 1 = 2Uρ. In contrast to homogeneous atomic gases considered in Refs. [10,11], where ρ is fixed and µ(ρ) should be calculated as an output parameter, in the triplon gas the chemical potential is fixed by the external magnetic field, while the density ρ = ρ(µ) should be calculated self consistently. In fact, in a system of triplons µ characterizes an additional direct contribution to the triplon energy due to the external field H and can be written as which can be interpreted as a chemical potential of the S z = −1 triplons.
The effects of disorder in dimerized quantum magnets in mean field approximations 6 The magnetization is proportional to the triplon density with ρ is defined from (8) as where µ 0 and µ 1 are given in (7) and the densities ρ 0 , ρ 1 must be calculated self consistently.
It is well known [24] that the disorder field leads to formation of a glassy fraction with the density ρ G . In this approximation each of ρ 1 and σ are presented as where ρ N and σ N are the normal and anomalous densities without disorder. In the YG method, based on HFB approximation, the following explicit relations can be obtained [18]: with the Bose distribution of Bogoliubov excitations f B ( For small momentum k the dispersion is linear, E k = ck, and the speed of the Bogoliubov mode The self energy ∆ is determined formally by the same equation as in the case when the disorder is neglected, The contribution from the disorder potential is hidden in the density of the glassy fraction where the double angle brackets mean the stochastic average. In general, the calculation of ρ G is rather complicated, however, for the δ -correlated disorder i.e. for the white noise, The effects of disorder in dimerized quantum magnets in mean field approximations 7 equation (18) is simplified as [10] The density of condensed fraction can be found by inserting (12) and (20) into the normalization condition (4). The result is In Eqs. (20) and (21) we introduced the dimensionless parameter R 0 as One can see from Eqs. (20) and (21) that the glassy fraction is proportional to the condensed one, The system of Eqs. (7), (8), (13)- (19) are the basic of YG approximation. An interesting quantity, crucial for determining the Bose glass phase, is the superfluid density, ρ s . In general it is defined as a partial density appearing as a response to a velocity boost whereP v is the total momentum of the system, dependent on the macroscopic velocity v. Referring the reader to original papers [10,11] we bring below analytical expression obtained there for ρ s in the case of white noise random potential Note that YG approach is valid for arbitrary strength of the interaction potential U, and for arbitrary strong disorder. For the weak interactions it leads to pioneering Huang-Meng approach [11], which will be extended to the "dirty triplons" in the next section.

Huang-Meng approximation
For completeness, we present here the results for the Huang-Meng approach, based on the so called Hartree Fock Popov (HFP) approximation which has been widely applied in the literature to describe the BEC of triplons [17,19]. The basic equations of this approach can be obtained by neglecting the anomalous density σ, which leads naturally to the single chemical potential µ = µ 0 = µ 1 . Namely, one finds from (7), (8) and (17) The effects of disorder in dimerized quantum magnets in mean field approximations 8 From these equations and (12) one obtains following main equations for the self energy ∆: where ρ N is formally given in (13), and ρ 0 is determined by the first equation in (27). The glassy fraction can be obtained from (20) in the linear approximation by R assuming weakness of interparticle interaction [10,11] where a s = Um/4π is the s -wave scattering length. Inserting (29) into (28) we can rewrite the former as To evaluate the densities one has to solve nonlinear algebraic equation (30), where ρ N is given formally by (13), with respect to ∆. Next, by inserting the result into (27) and (29) one obtains the density of condensed triplons ρ 0 and the glassy fraction ρ G , respectively. The total density can be evaluated then by the normalization condition ρ = ρ 0 + ρ N + ρ G . Equations (25), (26) for the superfluid density are formally the same in both approximations.
4. The shift of the critical temperature due to disorder and the T > T c regime It is well known that the critical temperature of BEC, T c for an ideal gas is given by: where ρ c is the total density of triplons near the critical temperature of BEC for pure system, with ζ(x) being the Riemann function. Equation (32) directly follows, from Eqs. (7), (8) or (28) by setting ρ N = ρ and ρ 0 = ρ G = 0. Clearly, any type of interaction is expected to modify T c . In general, these modifications are related to the interparticle interactions as well as to the disorder potential. Both approaches, considered here give a zero shift due to the boson-boson repulsion. However the shift due to the δ−correlated disorder (19), ∆T c = T c − T 0 c is given as [10,25] ∆T c T 0 where the dimensionless disorder parameter ν The effects of disorder in dimerized quantum magnets in mean field approximations 9 is introduced with the localization length For practical calculations we rewrite T c in Eq. (33), which is in a good agreement with perturbative estimates [2] as well as with Monte Carlo simulations [26], as an explicit function of effective mass m, the interaction strength U, critical magnetic field H c , disorder parameter ν, and external field H as follows: Now we proceed to consider the triplon density in the normal state in the T − T c ≫ ∆T c temperature range. The dirty bosons in the normal phase where the gauge symmetry is not broken, are yet poorly studied. For R = 0 with ρ 0 = ρ G = σ = 0 the triplon gas behaves like an "ideal gas" with an effective chemical potential µ eff , and the density [27] ρ Although µ eff is not accurately known it depends in general, on ρ, as well as on R. For the pure case MFA [27] gives µ eff (R = 0) = µ −2Uρ. The contribution from the disorder potential has been studied neither in YM nor in HM approaches. Therefore, to make the calculations self consistently, we have to use which yields the density ρ as a solution of the nonlinear equation (38).

Results and discussions
In the calculations below, the energies are measured in Kelvin, the mass in K −1 , the densities are dimensionless and the Bohr magneton is µ B = 0.671668 K/T. As to the strength of disorder potential R, it has units K −2 while the disorder parameter ν, defined in Eq. (34) is a number supposed to be less than one, ν < 1. As a material parameter, we use mean dimer-dimer distance in TlCuCl 3 r dd = 0.79 nm [19].
In Figure 1 we present as an example the total triplon density ρ(T ) for a clean and strongly disordered (ν = 0.45, see Eq.(34)) TlCuCl 3 , obtained in the YG approximation assuming that the total effect of the doping leads only to the randomness in the triplon subsystem.
The effects of disorder in dimerized quantum magnets in mean field approximations 10  The calculation of other quantities using the same assumption shows that the disorder leads to a decrease in the condensed and superfluid fractions, thereby increasing the glassy one. This tendency is quite natural, since the localization effects prevent particles from going into BEC. However, the increase in ρ G is so weak that along with ρ 0 the total number of triplons ρ is also decreased with increasing the strength of disorder potential R. Bearing in mind that ρ is proportional to the magnetization M, and ν is assumed to be approximately proportional to x, and comparing Fig.1 with the experimental data illustrated in Fig.2 one may conclude that the agreement between the theory and the experiment is unsatisfactory since the main features of the experimental results are not reproduced there. As it is seen in Fig.2 the disorder leads to an increase in the magnetization and, hence, in the total triplon density. This is accompanied by the decrease in the transition temperature. We therefore conclude that while the triplon gas can be considered similarly to atomic gases for which the considered meanfield approximations were developed, some further additional specific material -related properties of the dirty boson problem in quantum magnets must be taken into account.
First we note that the singlet -triplet excitation gap ∆ st , proportional to the critical field, H c , decreases under high pressure. This was experimentally observed in Ref. [28] for the pure spin system TlCuCl 3 . On the other hand it can be argued that the doping acts as a chemical pressure, which decreases H c . In fact, since the ionic radius of K + is smaller than that of Tl + , a partial substitution of Tl + ions with K + ions produces not only the exchange randomness, but also a compression of the crystal lattice. Thus the increase of the doping parameter, x, leads to decrease in H c which has indeed been observed experimentally [14,29,30]. Second, the disorder may increase the triplon effective mass The experimental low temperature magnetization in units of Bohr magneton per Cu ion of Tl 1−x K x CuCl 3 for various x in H = 7 T magnetic field obtained in Ref. [14].
thereby decreasing the critical temperature T c even when the gap decreases (similar effects were observed for helium in porous media [31,32]). Note that this effect manifests itself in different ways. For example, for the mixed compound IPACu(Cl x Br 1−x ) 3 the critical field, H c remains almost unchanged with varying x and then, abruptly becomes zero near the Cl-rich phase [33]. In another triplon-BEC compound, Ni(Cl 1−x Br x ) 2 -4SC(NH 2 ) 2 , it decreases by a factor of two when x changes from zero to 0.08 [30] although the physics of this decrease can be different from that in Tl 1−x K x CuCl 3 due to the fact that Br atomic radius is larger than the atomic radius of Cl. These effects of renormalization of the triplon spectrum by disorder can be considered similarly to the virtual crystal approximation in the simulations of disorder in solids, where the disorder is assumed to lead to a uniform change in the system parameters. The effects of disorder such as the appearance of the glassy phase with the density ρ G and related phenomena manifest themselves in addition to these uniform changes.
The phase diagram of Tl 1−x K x CuCl 3 in the (H, T ) plane was experimentally determined in Refs. [14,29] for various doping x, and the critical field H c was also estimated by extrapolation to zero temperature. In the present work the T c (H) dependence is given by Eq. (36). We made an attempt to least -square fit our parameters m and ν by using Eq. (36) to describe the experimental phase diagram. For simplicity we assume that interparticle interaction is not changed by doping, i.e. U = U(R = 0) = 313 K. The parameters obtained by this optimization are presented in Table 1.
Having fixed the input parameters for certain values of x, we are now in the position of recalculating the densities as well the magnetization to compare them with the Table 1. Optimized values of the input parameters of the model: the critical field, H c (taken from Ref. [14]), the disorder parameter ν, and the effective mass m for various doping x. The critical density, ρ c , the healing length, λ = 1/ √ 2mµ, the interparticle distance, d = 1/ρ 1/3 c and the localization length, L loc = d/ν are estimated at H = 7 T. It is assumed that the doping effects does not modify the Landé factor g and U .
x    Table 1. The total density of triplons is shown in Fig.3(d).
experiment. Figure 3 shows that the doping decreases ρ 0 and ρ s , and increases ρ G as it is expected due to the introduced disorder. Due to change of H c with x, the total density of triplons and hence the magnetization, now increases with increasing x in accordance with the experiment. One may conclude that the YG approach may well describe the effect of disorder to the magnetization, with the additional assumption of an x dependence of the effective mass and the critical field. One of the main characteristics of Bose condensed systems is the speed of the Bogoliubov mode c, defined here by the Eq. (16), which characterizes the propagation of collective excitations in the condensate. It is interesting to mention that the magnitude of c is large, being only an order of magnitude less than the speed of sound in the crystal. This is due to very small triplon mass in TlCuCl 3 . Clearly, disorder modifies the small-momentum excitation spectrum of the BEC. Estimates of such modification, ∆c = c−c 0 , where c 0 is the speed of the Bogoliubov mode for the pure system, that exist in the literature are controversial. For example, perturbative [1] and hydrodynamic [3] approaches give ∆c > 0, while ∆c < 0 was predicted in Refs. [4,5]. In Fig. 4 we present the corresponding speed for various doping parameters. It can be seen from comparison of Fig.4(a) and Fig.4(b) that both MFA approximations considered here show the decrease in c with increasing the disorder strength due to the localization effects. However, the effect of disorder is small leading to a less than 10 percent decrease in the speed of the Bogoliubov mode.
However, when the spectrum modification by disorder is also taken into account by a renormalization of the triplon mass and the gap, as close to the real situation, the dispersion of the sound-like mode in fixed magnetic field slightly increases with increasing disorder, reaches a maximum and then starts to decrease, (see Fig.4b). This behavior is caused by interplay between renormalization of the system parameters and localization effects. The former tends to increase c, e.g. by increasing µ, and therefore increasing the density, while the latter tends to decrease c, e.g. by decreasing the condensed fraction. Note that an increase in c with increasing the density was experimentally observed by Andrews et al. [34] for the BEC of sodium atoms. This interplay is illustrated in Fig.4(c) for ρ 0 and ρ s . It can be seen that uniform spectrum renormalization first leads to an antidepletion effect, increasing these quantities, while the localization effects impair the condensation and superfluidity.
We now consider the question about the existence of a pure Bose glass phase at T = 0, which, strictly speaking, should fulfill the following criteria [10,24,35]: (i) gapless in the excitation spectrum, (ii) insulating behavior, i.e. the superfluid fraction, ρ s = 0, (iii) finite compressibility, and (iv) finite density of states.
In 1970, Tachiki and Yamada [36] have shown that the Heisenberg-like Hamiltonian of s = 1/2 dimers can be rewritten as an effective bosonic Hamiltonian. Recently, Roscilde and Haas [37] generalized this bosonization procedure taking into account disorder and derived a Bose -Hubbard like Hamiltonian usually applied to study "dirty bosons" in optical lattices. Applying Fishers ideas [24] we may expect the formation of a pure Bose glass phase for doped magnets such as Tl 1−x K x CuCl 3 . Although Monte Carlo calculations [22,37] confirmed its existence the experimental confirmation is still a matter of debate [20,38]. We underline here that these Bose glass phases are localized out of the BEC phase, i.e. for H < H c . However, in the present work we have been mainly concentrating on the region with H ≥ H c where the gapless phase can be realized only within the BEC phase. For this case the definition of BG phase may be simplified  as a phase with ρ 0 = 0 and ρ s = 0 since the spectrum of the BEC is gapless by itself.
In searching for such a phase we studied ρ s and ρ 0 at T = 0 for various H ≥ H c and x and found no pure BG phase with ρ s = 0 as illustrated in Fig. 5. Note also that, as it is seen from   Table 1.
L loc > d. The physics of the possible BG phase for H < H c fields will be the subject of a separate study.

Conclusions
We reformulated and applied two existing mean field approximations for the "dirty boson" problem to study properties of Tl 1−x K x CuCl 3 quantum magnets. We showed that these approaches can qualitatively explain the magnetization data if a certain modification of the model parameters similar to the virtual crystal approximation usually applied to electron spectrum of alloys is taken into account. In fact, random bond effects in mixed dimerized magnetic compounds manifest themselves in a dual way: (i) by modification of internal parameters and (ii) by localization on random scatterers. Each of these effects could be studied separately in an appropriate theory, but they should be taken into account simultaneously for an adequate description of the measured magnetization data. Although the system becomes considerably disordered, the Bose-Einstein condensation does not support formation of the pure Bose glass phase. The random bonds lead to a nontrivial behavior of the sound-like mode speed: when H is fixed and x is experimentally varied, it increases for small x, reaches a maximum value and then decreases. While the speed of this mode was measured in dilute BEC of sodium atoms a long time ago [34], it has never been an intense focus of research in dimerized quantum magnets [40]. It could be systematically studied, for example, by measuring the dispersion relation of the Bogoliubov mode with inelastic neutron scattering techniques.

Introduction
The effects of disorder on the properties of Bose-Einstein condensates present interesting problems both for theoretical and experimental physics [1,2,3,4,5,6]. Disorder is important in various systems of real particles such as superfluid 4 He, cold atoms in optical lattices, and quasiparticles such as polaritons [7] and excitons [8]. These systems are well-suited for experimental studies, however the theory of disordered ensembles of interacting bosons is complex and there are essentially no exact solutions even in one dimension [9]. To approach this problem, Yukalov and Graham (YG) developed a self-consisted stochastic mean field approximation (MFA) [10] for Bose systems with arbitrarily strong interparticle repulsion and arbitrary strength of disorder potential. It was shown that, in general, the Bose system consists of following coexisting components: the condensate fraction, ρ 0 , the normal fraction ρ N , the glassy fraction ρ G , and, in addition, can be characterized by the superfluid density ρ s . In the limit of asymptotically weak interactions and disorder the known results, obtained in pioneering work by Huang and Meng [11] (HM) are reproduced by the YG theory. An interesting question here concerns the problem about the existence of a pure Bose glass (BG) phase, i.e. the phase where the condensate fraction is nonzero, while the superfluid fraction is not yet present. Note that Ref. [6] introduced an alternative definition of the gapless BG phase, having localized short-lived excitations and vanishing superfluid density with a continuous transition to the normal phase at finite temperature. Even without disorder, the condensate is depleted by particle-particle interactions and temperature. The inclusion of random fields depletes the condensate further and, possibly, creates the glassy fraction.
As it was understood recently, a new class of BEC can be provided by spinrelated quasiparticles in magnetic solids such as intensively pumped magnons [12] or triplons in the dimerized quantum magnets in the equilibrium [13]. In the magnets, the effect of disorder, which can be produced by admixing other chemical elements, can be rather strong to be seen in the physical properties such as the temperature-dependent magnetization. The so far most investigated compound showing BEC of triplons is TlCuCl 3 . To study the effect of disorder, solid solutions of quantum antiferromagnets TlCuCl 3 and KCuCl 3 , i.e. Tl 1−x K x CuCl 3 have been experimentally investigated recently [14,15,16] at low temperatures T . The zero-field ground states of TlCuCl 3 and KCuCl 3 are spin singlets with excitation gaps ∆ st = 7.1 K and ∆ st = 31.2 K, respectively and the magnetic excitations are spin triplets. Triplons arise in magnetic fields H > H c , where H c is defined by condition of closing the gap by the Zeeman splitting, that is ∆ st = gµ B H c , where g is the electron Landé factor and µ B is the Bohr magneton. In the mixture Tl 1−x K x CuCl 3 the induced magnetization M exhibits a cusplike minimum at a critical temperature T c (H) for fixed magnetic field H ≥ H c similarly to the parent compound, and can be successfully explained in terms of triplon BEC [17,18,19].
For a theoretical description it is natural to assume that for weak doping x ≪ 1 in the mixed system Tl 1−x K x CuCl 3 a small admixture of potassium forms a disorder potential. Consequently, the recently developed theories of "dirty bosons" [1,2,10,11] can be applied to study the BEC of triplons in Tl 1−x K x CuCl 3 . Here the following natural questions arise. For example, what is the correspondence between admixing parameter x and the properties of the disorder potential ? What are the experimental consequences of the disorder ? Yamada et al. [15] analyzed the electron spin resonance spectrum in Tl 1−x K x CuCl 3 and concluded that there is a Bose glass -BEC transition near a critical magnetic field. Although this interpretation might need a further analysis (see discussion in Refs. [20]) it would be interesting to study the influence of the glassy phase, or more exactly, of the glassy fraction ρ G on the magnetization. Note that even the existence of a pure Bose glass phase still is a matter of debate even in theoretical approaches. For example, it may be predicted by the approach used by Huang and Meng [11] if one extends their formulas from weak disorder to a strong one. On the other hand, no pure Bose glass was found in Monte -Carlo simulations [21] for atomic gases, but predicted for triplons at T = 0 by Nohadani et al. [22].
Here we develop a theory of the disorder effects on the BEC of triplons taking Tl 1−x K x CuCl 3 as a prototype for studies of specific properties. For example, in atomic gases considered in Refs. [10,11], the chemical potential µ is determined self-consistently with fixed number of atoms, while in the triplon gas the chemical potential is a given external parameter controlled by the applied magnetic field and the number of triplons is conserved in the thermodynamic limit. To clarify the terms, we underline that the number of magnons may vary but that of triplons may be tuned and kept fixed, which makes possible the BEC of the latter.
The paper is organized as follows. In Sections II and III we outline the YG and HM approaches valid only for T ≤ T c and extend it for the triplon system. The shift of T c due to disorder and the normal phase properties will be discussed in Section IV. Our numerical results will be presented in Section V. Conclusions will summarize the results of this work.

Yukalov-Graham approximation for disordered triplons
In the following we reformulate the Yukalov-Graham approximation to the triplon system with arbitrary disorder. The Hamiltonian operator of triplons with contact interaction and implemented disorder potential V (r) is given by where ψ(r) is the bosonic field operator, U is the interparticle interaction strength, and K is the kinetic energy operator which defines the bare triplon dispersion ε k . Since the triplon BEC occurs in solids, we perform integration over the unit cell of the crystal with the corresponding momenta defined in the first Brillouin zone. Below we assume that the bare spectrum remains coherent in the presence of disorder and consider it as a simple isotropic one: ε k = k 2 /2m, where m is the triplon effective mass. The distribution of random fields is assumed to be zero -centered, V (r) = 0, and the correlation function Here and below we adopt the units k B ≡ 1, ≡ 1, and V ≡ 1 if not stated otherwise for the unit cell volume.
To describe Bose condensed system where the global gauge symmetry is broken, one employs the Bogoliubov shift: where the condensate density ρ 0 (r) is constant for the homogeneous system, ρ 0 (r) ≡ ρ 0 . Since by the definition the average of ψ † (r)ψ(r) is the total number of particles: with the density of triplons per unit cell ρ = N/V , from the normalization condition one immediately obtains Therefore the field operator ψ 1 (r) determines the density of uncondensed particles. The YG approximation is formulated in representative ensemble formalism, which includes two Lagrange multipliers, µ 0 and µ 1 , defined as: where Ω is the grand thermodynamic potential. It was shown that disorder would not change the explicit expressions for chemical potentials, obtained earlier [23] in Hartree-Fock-Bogoliubov (HFB) approximation without disorder, where σ = 1 V V d 3 r ψ 1 (r)ψ 1 (r) is the anomalous density. The total system chemical potential mu, related to the total number of particles as N = −∂Ω/∂µ, is determined by Clearly, when the gauge symmetry is not broken, i.e. ρ 0 = 0, σ = 0, ρ 1 = ρ, both µ 0 and µ 1 coincide giving µ = µ 1 = 2Uρ. In contrast to homogeneous atomic gases considered in Refs. [10,11], where ρ is fixed and µ(ρ) should be calculated as an output parameter, in the triplon gas the chemical potential is fixed by the external magnetic field, while the density ρ = ρ(µ) should be calculated self consistently. In fact, in a system of triplons µ characterizes an additional direct contribution to the triplon energy due to the external field H and can be written as which can be interpreted as a chemical potential of the S z = −1 triplons.
The effects of disorder in dimerized quantum magnets in mean field approximations 6 The magnetization is proportional to the triplon density with ρ is defined from (8) as where µ 0 and µ 1 are given in (7) and the densities ρ 0 , ρ 1 must be calculated self consistently.
It is well known [24] that the disorder field leads to formation of a glassy fraction with the density ρ G . In this approximation each of ρ 1 and σ are presented as where ρ N and σ N are the normal and anomalous densities without disorder. In the YG method, based on HFB approximation, the following explicit relations can be obtained [18]: with the Bose distribution of Bogoliubov excitations f B ( For small momentum k the dispersion is linear, E k = ck, and the speed of the Bogoliubov mode The self energy ∆ is determined formally by the same equation as in the case when the disorder is neglected, The contribution from the disorder potential is hidden in the density of the glassy fraction where the double angle brackets mean the stochastic average. In general, the calculation of ρ G is rather complicated, however, for the δ -correlated disorder i.e. for the white noise, The effects of disorder in dimerized quantum magnets in mean field approximations 7 equation (18) is simplified as [10] The density of condensed fraction can be found by inserting (12) and (20) into the normalization condition (4). The result is In Eqs. (20) and (21) we introduced the dimensionless parameter R 0 as One can see from Eqs. (20) and (21) that the glassy fraction is proportional to the condensed one, The system of Eqs. (7), (8), (13)- (19) are the basic of YG approximation. An interesting quantity, crucial for determining the Bose glass phase, is the superfluid density, ρ s . In general it is defined as a partial density appearing as a response to a velocity boost whereP v is the total momentum of the system, dependent on the macroscopic velocity v. Referring the reader to original papers [10,11] we bring below analytical expression obtained there for ρ s in the case of white noise random potential Note that YG approach is valid for arbitrary strength of the interaction potential U, and for arbitrary strong disorder. For the weak interactions it leads to pioneering Huang-Meng approach [11], which will be extended to the "dirty triplons" in the next section.

Huang-Meng approximation
For completeness, we present here the results for the Huang-Meng approach, based on the so called Hartree Fock Popov (HFP) approximation which has been widely applied in the literature to describe the BEC of triplons [17,19]. The basic equations of this approach can be obtained by neglecting the anomalous density σ, which leads naturally to the single chemical potential µ = µ 0 = µ 1 . Namely, one finds from (7), (8) and (17) The effects of disorder in dimerized quantum magnets in mean field approximations 8 From these equations and (12) one obtains following main equations for the self energy ∆: where ρ N is formally given in (13), and ρ 0 is determined by the first equation in (27). The glassy fraction can be obtained from (20) in the linear approximation by R assuming weakness of interparticle interaction [10,11] where a s = Um/4π is the s -wave scattering length. Inserting (29) into (28) we can rewrite the former as To evaluate the densities one has to solve nonlinear algebraic equation (30), where ρ N is given formally by (13), with respect to ∆. Next, by inserting the result into (27) and (29) one obtains the density of condensed triplons ρ 0 and the glassy fraction ρ G , respectively. The total density can be evaluated then by the normalization condition ρ = ρ 0 + ρ N + ρ G . Equations (25), (26) for the superfluid density are formally the same in both approximations.

The shift of the critical temperature due to disorder and the T > T c regime
It is well known that the critical temperature of BEC, T c for an ideal gas is given by: where ρ c is the total density of triplons near the critical temperature of BEC for pure system, with ζ(x) being the Riemann function. Equation (32) directly follows, from Eqs. (7), (8) or (28) by setting ρ N = ρ and ρ 0 = ρ G = 0. Clearly, any type of interaction is expected to modify T c . In general, these modifications are related to the interparticle interactions as well as to the disorder potential. Both approaches, considered here give a zero shift due to the boson-boson repulsion. However the shift due to the δ−correlated disorder (19), ∆T c = T c − T 0 c is given as [10,25] ∆T c T 0 where the dimensionless disorder parameter ν The effects of disorder in dimerized quantum magnets in mean field approximations 9 is introduced with the localization length For practical calculations we rewrite T c in Eq. (33), which is in a good agreement with perturbative estimates [2] as well as with Monte Carlo simulations [26], as an explicit function of effective mass m, the interaction strength U, critical magnetic field H c , disorder parameter ν, and external field H as follows: Now we proceed to consider the triplon density in the normal state in the T − T c ≫ ∆T c temperature range. The dirty bosons in the normal phase where the gauge symmetry is not broken, are yet poorly studied. For R = 0 with ρ 0 = ρ G = σ = 0 the triplon gas behaves like an "ideal gas" with an effective chemical potential µ eff , and the density [27] ρ Although µ eff is not accurately known it depends in general, on ρ, as well as on R. For the pure case MFA [27] gives µ eff (R = 0) = µ −2Uρ. The contribution from the disorder potential has been studied neither in YM nor in HM approaches. Therefore, to make the calculations self consistently, we have to use which yields the density ρ as a solution of the nonlinear equation (38).

Results and discussions
In the calculations below, the energies are measured in Kelvin, the mass in K −1 , the densities are dimensionless and the Bohr magneton is µ B = 0.671668 K/T. As to the strength of disorder potential R, it has units K −2 while the disorder parameter ν, defined in Eq. (34) is a number supposed to be less than one, ν < 1. As a material parameter, we use mean dimer-dimer distance in TlCuCl 3 r dd = 0.79 nm [19].
In Figure 1 we present as an example the total triplon density ρ(T ) for a clean and strongly disordered (ν = 0.45, see Eq.(34)) TlCuCl 3 , obtained in the YG approximation assuming that the total effect of the doping leads only to the randomness in the triplon subsystem.  The calculation of other quantities using the same assumption shows that the disorder leads to a decrease in the condensed and superfluid fractions, thereby increasing the glassy one. This tendency is quite natural, since the localization effects prevent particles from going into BEC. However, the increase in ρ G is so weak that along with ρ 0 the total number of triplons ρ is also decreased with increasing the strength of disorder potential R. Bearing in mind that ρ is proportional to the magnetization M, and ν is assumed to be approximately proportional to x, and comparing Fig.1 with the experimental data illustrated in Fig.2 one may conclude that the agreement between the theory and the experiment is unsatisfactory since the main features of the experimental results are not reproduced there. As it is seen in Fig.2 the disorder leads to an increase in the magnetization and, hence, in the total triplon density. This is accompanied by the decrease in the transition temperature. We therefore conclude that while the triplon gas can be considered similarly to atomic gases for which the considered meanfield approximations were developed, some further additional specific material -related properties of the dirty boson problem in quantum magnets must be taken into account.
First we note that the singlet -triplet excitation gap ∆ st , proportional to the critical field, H c , decreases under high pressure. This was experimentally observed in Ref. [28] for the pure spin system TlCuCl 3 . On the other hand it can be argued that the doping acts as a chemical pressure, which decreases H c . In fact, since the ionic radius of K + is smaller than that of Tl + , a partial substitution of Tl + ions with K + ions produces not only the exchange randomness, but also a compression of the crystal lattice. Thus the increase of the doping parameter, x, leads to decrease in H c which has indeed been observed experimentally [14,29,30]. Second, the disorder may increase the triplon effective mass The experimental low temperature magnetization in units of Bohr magneton per Cu ion of Tl 1−x K x CuCl 3 for various x in H = 7 T magnetic field obtained in Ref. [14].
thereby decreasing the critical temperature T c even when the gap decreases (similar effects were observed for helium in porous media [31,32]). Note that this effect manifests itself in different ways. For example, for the mixed compound IPACu(Cl x Br 1−x ) 3 the critical field, H c remains almost unchanged with varying x and then, abruptly becomes zero near the Cl-rich phase [33]. In another triplon-BEC compound, Ni(Cl 1−x Br x ) 2 -4SC(NH 2 ) 2 , it decreases by a factor of two when x changes from zero to 0.08 [30] although the physics of this decrease can be different from that in Tl 1−x K x CuCl 3 due to the fact that Br atomic radius is larger than the atomic radius of Cl. These effects of renormalization of the triplon spectrum by disorder can be considered similarly to the virtual crystal approximation in the simulations of disorder in solids, where the disorder is assumed to lead to a uniform change in the system parameters. The effects of disorder such as the appearance of the glassy phase with the density ρ G and related phenomena manifest themselves in addition to these uniform changes.
The phase diagram of Tl 1−x K x CuCl 3 in the (H, T ) plane was experimentally determined in Refs. [14,29] for various doping x, and the critical field H c was also estimated by extrapolation to zero temperature. In the present work the T c (H) dependence is given by Eq. (36). We made an attempt to least -square fit our parameters m and ν by using Eq. (36) to describe the experimental phase diagram. For simplicity we assume that interparticle interaction is not changed by doping, i.e. U = U(R = 0) = 313 K. The parameters obtained by this optimization are presented in Table 1.
Having fixed the input parameters for certain values of x, we are now in the position of recalculating the densities as well the magnetization to compare them with the Table 1. Optimized values of the input parameters of the model: the critical field, H c (taken from Ref. [14]), the disorder parameter ν, and the effective mass m for various doping x. The critical density, ρ c , the healing length, λ = 1/ √ 2mµ, the interparticle distance, d = 1/ρ 1/3 c and the localization length, L loc = d/ν are estimated at H = 7 T. It is assumed that the doping effects does not modify the Landé factor g and U .
x    Table 1. The total density of triplons is shown in Fig.3(d).
experiment. Figure 3 shows that the doping decreases ρ 0 and ρ s , and increases ρ G as it is expected due to the introduced disorder. Due to change of H c with x, the total density of triplons and hence the magnetization, now increases with increasing x in accordance with the experiment. One may conclude that the YG approach may well describe the effect of disorder to the magnetization, with the additional assumption of an x dependence of the effective mass and the critical field. One of the main characteristics of Bose condensed systems is the speed of the Bogoliubov mode c, defined here by the Eq. (16), which characterizes the propagation of collective excitations in the condensate. It is interesting to mention that the magnitude of c is large, being only an order of magnitude less than the speed of sound in the crystal. This is due to very small triplon mass in TlCuCl 3 . Clearly, disorder modifies the small-momentum excitation spectrum of the BEC. Estimates of such modification, ∆c = c−c 0 , where c 0 is the speed of the Bogoliubov mode for the pure system, that exist in the literature are controversial. For example, perturbative [1] and hydrodynamic [3] approaches give ∆c > 0, while ∆c < 0 was predicted in Refs. [4,5]. In Fig. 4 we present the corresponding speed for various doping parameters. It can be seen from comparison of Fig.4(a) and Fig.4(b) that both MFA approximations considered here show the decrease in c with increasing the disorder strength due to the localization effects. However, the effect of disorder is small leading to a less than 10 percent decrease in the speed of the Bogoliubov mode.
However, when the spectrum modification by disorder is also taken into account by a renormalization of the triplon mass and the gap, as close to the real situation, the dispersion of the sound-like mode in fixed magnetic field slightly increases with increasing disorder, reaches a maximum and then starts to decrease, (see Fig.4b). This behavior is caused by interplay between renormalization of the system parameters and localization effects. The former tends to increase c, e.g. by increasing µ, and therefore increasing the density, while the latter tends to decrease c, e.g. by decreasing the condensed fraction. Note that an increase in c with increasing the density was experimentally observed by Andrews et al. [34] for the BEC of sodium atoms. This interplay is illustrated in Fig.4(c) for ρ 0 and ρ s . It can be seen that uniform spectrum renormalization first leads to an antidepletion effect, increasing these quantities, while the localization effects impair the condensation and superfluidity.
We now consider the question about the existence of a pure Bose glass phase at T = 0, which, strictly speaking, should fulfill the following criteria [10,24,35]: (i) gapless in the excitation spectrum, (ii) insulating behavior, i.e. the superfluid fraction, ρ s = 0, (iii) finite compressibility, and (iv) finite density of states.
In 1970, Tachiki and Yamada [36] have shown that the Heisenberg-like Hamiltonian of s = 1/2 dimers can be rewritten as an effective bosonic Hamiltonian. Recently, Roscilde and Haas [37] generalized this bosonization procedure taking into account disorder and derived a Bose -Hubbard like Hamiltonian usually applied to study "dirty bosons" in optical lattices. Applying Fishers ideas [24] we may expect the formation of a pure Bose glass phase for doped magnets such as Tl 1−x K x CuCl 3 . Although Monte Carlo calculations [22,37] confirmed its existence the experimental confirmation is still a matter of debate [20,38]. We underline here that these Bose glass phases are localized out of the BEC phase, i.e. for H < H c . However, in the present work we have been mainly concentrating on the region with H ≥ H c where the gapless phase can be realized only within the BEC phase. For this case the definition of BG phase may be simplified  T as a function of doping parameter x, taking into account solely renormalization of the system parameters in Table 1 in the YG (solid line) and the HM (dashed line) approaches. (b) The same as in Fig. 4(a), now with effects of disorder taken into account. (c) The superfluid and condensed fractions (as marked near the plots) in the YG (solid lines) and the HM (dashed lines) approximations with the renormalized bare spectrum parameters.
as a phase with ρ 0 = 0 and ρ s = 0 since the spectrum of the BEC is gapless by itself.
In searching for such a phase we studied ρ s and ρ 0 at T = 0 for various H ≥ H c and x and found no pure BG phase with ρ s = 0 as illustrated in Fig. 5. Note also that, as it is seen from Table 1, for moderate values of x considered here, the localization length, i.e the mean free path [39] is larger than interparticle distance, Upper panel corresponds to the YG approximation and lower panel corresponds to the HM approximation presented here for comparison. Graphs in plots (a), (c) were calculated without bare spectrum renormalization, while graphs in plots (b), (d) take into account spectrum renormalization presented in Table 1.
L loc > d. The physics of the possible BG phase for H < H c fields will be the subject of a separate study.

Conclusions
We reformulated and applied two existing mean field approximations for the "dirty boson" problem to study properties of Tl 1−x K x CuCl 3 quantum magnets. We showed that these approaches can qualitatively explain the magnetization data if a certain modification of the model parameters similar to the virtual crystal approximation usually applied to electron spectrum of alloys is taken into account. In fact, random bond effects in mixed dimerized magnetic compounds manifest themselves in a dual way: (i) by modification of internal parameters and (ii) by localization on random scatterers. Each of these effects could be studied separately in an appropriate theory, but they should be taken into account simultaneously for an adequate description of the measured magnetization data. Although the system becomes considerably disordered, the Bose-Einstein condensation does not support formation of the pure Bose glass phase. The random bonds lead to a nontrivial behavior of the sound-like mode speed: when H is fixed and x is experimentally varied, it increases for small x, reaches a maximum value and then decreases. While the speed of this mode was measured in dilute BEC of sodium atoms a long time ago [34], it has never been an intense focus of research in dimerized quantum magnets [40]. It could be systematically studied, for example, by measuring the dispersion relation of the Bogoliubov mode with inelastic neutron scattering techniques.