Random-field-induced order in bosonic t-J model

In the present paper, we shall study effect of a random quenched external field for spin order and also multiple Bose-Einstein condensation (BEC). This system is realized by the cold atomic gases in an optical lattice. In particular, we are interested in the strong-repulsion region of two-component gases for which the bosonic t-J model is a good effective model. In the bosonic t-J model, a long-range order of the pseudo-spin and also BEC of atoms appear quite naturally as in the fermion t-J model for the high-temperature superconducting materials. Random Raman scattering between two internal states of a single atom plays a role of the random external field, and we study its effects on the pseudo-spin order and the BEC by means of quantum Monte-Carlo simulations. The random external field breaks a continuous U(1) symmetry existing in the original bosonic t-J model and it induces new orders named random-field-induced order (RFIO). We show a phase diagram of the bosonic t-J model with the random external magnetic field and study the robustness of the RFIO states. We also study topological excitations like vortices and domain wall in the RFIO state. Finally we point out the possibility of a quantum bit by the RFIO.


Introduction
Quenched disorder plays very important role in condensed matter physics. A prominent example is the Anderson localization that predicts all quantum states are localized in one and two spatial dimensions if interactions between particles can be neglected [1]. It has been proved rigorously that quenched disorder destroys ordered states and rounds singularity of phase transitions [2]. However recently, a counter-intuitive possibility was pointed out and examined. That is, a quenched disorder generates a new ordered state, which is different from the original one, if that quenched disorder breaks a continuous symmetry of the original system without the quenched disorder. This possibility was first studied in a classical XY spin model coupled with a random external magnetic field [3] and then the resultant order is called random-field-induced order (RFIO). Shortly after the proposal, it was shown that such phenomena of the RFIO can be observed by experiments on systems of ultra-cold bosonic atoms of multiple-internal states [4]. In Bose atomic gases, a Bose-Einstein condensed (BEC) state is the genuine ordered state, and a random Raman scattering of the internal states of the atom plays a role of a quenched random external magnetic field. Then, it is expected that the RFIO can be observed in the ultra-cold atomic gas systems. There are other interesting works on the RFIO [5].
In this paper, we shall investigate the RFIO in detail for the two-component cold bosonic gases in a square optical lattice (OL). In particular, we consider the bosonic t-J model, which is a low-energy effective model for the Bose-Hubbard model in the strong-repulsion limit [6,7]. This model exhibits both the pseudo-spin order and the BECs and is therefore suitable for study of the RFIO. We employ quantum Monte-Carlo (MC) simulations that take into account all of fluctuations. Thus the present study is in sharp contrast to the previous ones that used estimation of classical energy of the XY spin configurations for observing a possible RFIO [3] and a mean-field theory with Gross-Pitaevskii equations for the BEC of the RFIO [4]. Furthermore we will investigate behavior of low-energy topological excitations such as vortices and domain walls in the RFIO states, and reveal interesting properties of them. Finite-temperature (T ) phase diagram is also obtained, which is useful for discussion on the robustness of the RFIO states.
The present paper is organized as follows. In Sec.2, we introduce and explain the bosonic t-J model with a random external field. Path-integral quantization using the slave-particle representation is explained. Effective field theory for the pseudo-spin and BECs is derived by integrating out amplitude degrees of freedom of the slave-particle field variables. In Sec.3, a replica mean-field theory is applied to the effective field theory and effects of the random field are studied. This study clearly shows how the RFIO appears as a result of the random field with a moderate fluctuation. In Sec.4, the results of the numerical simulations are given. Phase diagrams at vanishing T as well as finite-T are obtained. Various correlation functions, which are used for the identification of the orders are shown. In Sec.5., topological excitations like vortices and domain wall are investigated numerically. Properties of these excitations are discussed from the view point of the RFIO. Section 6 is devoted for conclusion.

Bosonic t-J model with random Rabi coupling
The bosonic t-J model, which describes dynamics of two internal states of a boson, which we call a-and b-boson for simplicity, in a square OL, is defined by the following Hamiltonian [7], where a † i (a i ) and b † i (b i ) are boson creation (destruction) operators at site i of the square lattice and t a and t b are the hopping amplitude between the nearest-neighbor (NN) sites. Pseudo-spin operator S i is given as and σ is the Pauli spin matrix. In the t-J model, the doubly-occupied state is excluded at each site and density of atoms is controlled by the chemical potential µ. It was probed that the Hamiltonian H tJ is derived from the Bose-Hubbard model in the strong one-site repulsion limit by integrating out multiple-particle states, and the exchange couplings J xy and J z are related with the intra and inter-repulsions between atoms. In the present study, however, we shall treat these parameters as free ones because the system H tJ might be derived from a Bose-gas system on the Lieb lattice that is a bosonic counterpart of the d-p model for the strongly-correlated electron systems.
The term H V in Eq.(3) controls density fluctuations of atoms at each site from the mean valueρ ai andρ bi . This term is expected to appear naturally for describing practical phenomena in experiments at low energies and therefore we explicitly added it to the Hamiltonian.
In the later discussion, we shall mostly consider the case J z = 0, which corresponds to the case of the equal intra and inter-species repulsions in the Bose-Hubbard model. The bosonic t-J model without the random external field was studied by both the numerical and analytical methods in the previous papers and its phase diagram has been clarified [7,6]. For sufficiently large J xy , a ferromagnetic (FM) state of the pseudospin appears, whereas for sufficiently large hopping amplitude t a and t b , BECs of the atoms form. As the anti-ferromagnetic coupling J z is increased, supersolid forms for a small but finite parameter region. For sufficiently large J z compared with t a , t b and J xy , a solid state with the checkerboard density pattern appears as the lowest-energy state.
In the present study, we add the following terms that describe the quenched random external fields, where J x i and J y i take random real variables with the vanishing mean value. For the practical numerical study, we use the following distribution function P (J i ), where σ x(y) are positive parameters. In the cold atomic systems, the above terms are realized by the Rabi oscillation with the Raman laser of a random complex amplitude . It is expected that the complex Raman amplitude Ω is realized experimentally using speckle laser light [8].
It should be remarked that the U(1)×U(1) symmetry of H EtJ in Eq.(1), i.e., (a i , b i ) → (e iα a i , e iβ b i ) with arbitrary constants α and β, is preserved only for the case σ x = σ y , otherwise the quenched disorder, J x i and J y i , explicitly breaks the symmetry as U(1)×U(1)→U(1)×Z 2 .

Numerical methods:Path-integral Monte-Carlo simulations
In order to study the model H T by means of quantum MC simulations, we use the path-integral method with the slave particle description for the local constraint of the t-J model. The boson creation operators are expressed by the slave particle operators φ ai , φ bi and φ hi as follows, and physical state of the slave particle |P hys must satisfy Then the partition function for the system H T is given by where τ is the imaginary time,φ αi = dφ αi dτ and H T is expressed in terms of the slave particles by using Eq. (6). For the path integral in Eq.(8), the local constraint φ ai φ ai +φ bi φ bi +φ hi φ hi = 1 can be imposed by using a Lagrange multiplier field λ i (τ ), To obtain a positive-definite action for carrying out the path-integral MC simulation, we parameterize the fields as φ αi = √ ρ αi e iω αi (α = a, b, h), and analytically calculate the integral over the amplitudes ρ αi . By the term H V in H T , the integration can be carried out in powers of the density fluctuations, δρ αi = ρ αi −ρ αi . As a result, the Berry phase αφαiφαi generates terms like 1 V 0 α (ω αi + λ i ) 2 in the action. For the practical numerical calculation, we introduce a lattice for the imaginary time direction, and we denote a site of three-dimensional (3D) cubic space-time lattice r. With the resultant action on the lattice A Lxy , the partition function is given by with where and In Eqs.(11) ∼ (14), dynamical variables are θ sr = ω ar − ω br , θ ar = ω ar − ω hr , θ br = ω br − ω hr , and parameters are related to the original ones as, where ∆τ is the lattice spacing of the imaginary time. Note that c τ , · · · ,J y are all dimensionless. (We have puth = 1.) Please notice that the quenched disorder variables J x(y) i are independent of the imaginary time τ . There are comments on the derivation of A Lxy and advantages of the MC simulation on it. On performing the path integral of ρ αi , the higher-order terms of the fluctuations δρ αi are ignored, e.g., in the hopping term, The above approximation is legitimate for δρ/ρ ≪ 1 as in the experiments of largē ρ [9] or small δρ. In the previous paper [10], we studied the non-random case rather in detail and verified that results obtained by the MC simulations on A Lxy are in good agreement with those obtained by the Gross-Pitaevskii theory. Furthermore in the previous paper [11], we studied the phase diagram of H EtJ with a finite J z by using A Lxy and found that the supersolid state forms in certain parameter region. The obtained phase diagram is in agreement with that of the two-component Bose-Hubbard model, which was obtained by using the MC simulation with the worm algorithm [12], although the case of commensurate filling factors was studied there. One advantage of the present MC method for studying the bosonic t-J model is its rapid convergence, and therefore the large-scale MC simulation is possible. Furthermore, various correlation functions as well as the density of topological excitations can be calculated accurately, as we show in sections 4 and 5.

Replica mean-field theory
Before going into the numerical calculations, we briefly study the model given by Eq.(10) by means of the replica methods. In particular, we are interested in the case of the single-component random field likeJ x = 0 andJ y = 0, and see how the order of S y shows up whereas that of S x does not. For simplicity, we shall consider the case of the total filling factor =1, i.e., the filling factor of each particle is 1/2, and focus on the pseudo-spin symmetry, though the extension to the case with a finite hole density is rather straightforward.
In the replica method studying effects of quenched random variables, a replica index ν = 1, 2, · · · , n is introduced for each dynamical variable. In the present system, ω αi → ω ν αi (α = a, b, h), and the partition function of the replica system [Z n ] is given by, whereω ν ai = dω ν ai dτ , etc, and [· · ·] denotes average over the random variablesJ In [Z n ] in Eq.(17), the integration overJ x i can be carried out readily to obtain, The nonlocal terms in Eq.(18) can be reduced to local ones by using a Hubbard-Storatonovich transformation with auxiliary fields m i (τ ) as We also apply a mean-field theory (MFT) for the spin part of A ν S in Eq.(17) as In this MFT, the partition function of the replica system [Z n ] MFT is given as In Eq.(21), the integration of ω ν αi can be carried out to obtain a Ginzburg-Landau theory (GL theory) for the pseudo-spin order. To this end, we use the following on-site Green functions as we consider the system at sufficiently low temperature, Then where γ ν i = m i + 4C 1 S xν . We are interested in the replica-symmetric solution and set ν S xν = n S x , etc. We also introduce a cutoff β = 1/(k B T ) for the integral of the imaginary time τ . Then the integration over m i can be done to obtain Finally we obtain the effective potential, V Rep , by taking the limit n → 0, It is obvious that two limits, β → ∞ and n → 0, are not interchangeable in V Rep given by Eq.(26). For a finite σ x and at finite temperature, the last term in V Rep (26) vanishes for n → 0 and both S x and S y can have a nonvanishing value for C 1 > V 0 2 , i.e., in the case in which the spin interaction J dominates the suppression of the density fluctuations V 0 . On the other hand for the case of large fluctuation of the random field and then S x does not condense whereas S y does. Similarly for T → 0, the last term of V Rep in Eq.(26) gives a finite contribution for any nonvanishing σ x , and S x does not condense. The above results seem interesting but they are obtained by the MFT. For example, the assumption of the constant mean field S x is not correct for σ x → ∞. Therefore more reliable studies are welcome. The numerical calculations in the subsequent sections give reliable results and reveal detailed properties of the RFIO.

Phase diagrams at low temperature
In this and subsequent sections, we shall show the results obtained by means of the numerical MC simulations. Model is defined by Eqs.(10)∼ (14), and the local-update MC simulation was used for calculation of physical quantities for fixed random variables The standard Metropolis algorithm [13] was used for the local update. For the local update of the angle variables θ αi , random variables ∆θ used for generating a candidate of a new variable θ new = θ old + ∆θ was chosen in the range |∆θ| < π 3 . Typical sweeps for the thermalization is 100000 and for the measurement is (20000)×(10 samples). Typical acceptance ratio is 40% ∼ 50%, and errors were estimated from 10 samples by the jackknife methods [14]. We first show the phase diagram of the system without the random field, which was obtained in the previous study [7]. For the case of t a = t b and at T = 0 ‡, there are three phases, phase with no long-range order (LRO), FM phase and phase of double BECs, which we often denote 2SF, accompanying the FM order. See Fig.1. In particular in the states with the FM order, the pseudo-spin (S x , S y ) has a LRO in an arbitrary direction, i.e., S y / S x = tan θ s with an arbitrary angle θ s . This is observed by calculating the correlation functions G x S (r) and G y S (r) defined by G x(y) where L is the linear size of the 3D lattice, and sites r 0 and r 0 +r are located in the same spatial 2D lattice, i.e., G x(y) S (r) is an equal-time correlator. We put the lattice size of the imaginary-time direction N τ = the lattice size of the spatial direction L. The angle θ takes various values depending on initial configurations and random variables used local updates in the MC simulations. This result comes from the U(1) symmetry of the pseudo-spin rotation in the system of the action A Lτ + A L .
In the present paper, we shall study the system in the random fields A q in Eq. (14). We first consider the caseJ y i = 0. In this single-component external system, the U(1) spin symmetry is reduced to the Z 2 symmetry of the Ising type (S x , S y ) → −(S x , S y ). This fact implies that there might exist a preferred direction in the pseudo-spin order. This is actually the case as we show shortly.
In the numerical studies, we first determine the random variables {J x i } according to the distribution P (J x ) in Eq. (5). To this end, we used the box-Muller methods. Then the MC simulation is carried out for the system with the fixed {J x i } by the local update of ω αr and λ r . Final physical quantities are obtained by averaging calculated quantities ‡ More precisely in the MC simulations, temperature of the system T is given as kBT V0 = cτ Nτ , where N τ is the lattice size in the imaginary-time direction. Then the system at T → 0 is realized as N τ → ∞.  for each sample over 5 ∼ 10 {J x i } samples. Phase boundary is determined by calculating the internal energy E and the specific heat C, which are defined as where the mean value · includes the average over samples of the random fields as we explained above. In order to identify physical properties of each phase, we also calculate the boson correlation functions besides the pseudo-spin one in Eq.(28), r 0 e iθar 0 e −iθ a,r 0 +r , where, as in G x(y) S (r), sites r 0 and r 0 + r are located in the same spatial 2D lattice. In Fig.1, we exhibit the obtained phase diagram for σ x = 0.3 and σ y = 0, which corresponds toJ y i = 0. Typical behaviors of E and C of various system sizes, which are used for identification of the phase boundaries, are also shown in Fig.2. The calculations indicate that all phase transitions are of second order. The locations of the phase boundaries are almost the same with the ones of the original bosonic t-J model. However the FM order is replaced by the RFIO, as the pseudo-spin correlation function in Fig.3 indicates that only the y-component of the pseudo-spin has a LRO and the correlation of the x-component vanishes quite rapidly as a function of r. This result means that the relative phase of the condensations of the a-and b-boson operators has definite values θ s ≃ ± π 2 . To verify this, we measured the relative phase at each site and the result is shown in Fig.4. The magnetization in the phase of RFIO is slightly smaller than that in the original t-J model. The above observation is in good agreement with the previous studies of the related systems such as the classical XY spin model and two component BEC of the cold atoms [3,4].  It is interesting to see how the phase diagram is changed when both components of the random field are turned on. It is not so difficult to show that the U(1) symmetry is restored for the case σ x = σ y , and then the relative phase θ s = θ a − θ b takes an arbitrary value. From this observation, one may expect that the range of θ is expanded for nonvanishing σ x and σ y depending on the ratio σ x /σ y . However this is not the case. We studied the cases with various value of σ x /σ y , and found that θ s takes ± π 2 (0 or π) for σ x /σ y > 1 (< 1). For a typical example, see Fig.5.
Let us perform the "Gedanken experiment" in which the parameter V 0 is varied with the other parameters fixed. For larger V 0 , fluctuations of the densities of atoms in each site are suppressed and then fluctuations of the phase degrees of freedom of the boson operators are enhanced. In fact from the action in Eqs. (12) and (15), it is seen that the phases ω αr vary rapidly in the τ -direction for small c τ , even though their order are generated in the spatial direction for a sufficiently large hopping amplitude and a spin-exchange coupling. We show the results of the numerical study in Fig.6. Phase diagram is given in the ( 1 V 0 − σ x ) plane for C 1 = 2.0 and C a 3 = C b 3 = 0.2. In the 3D region of smaller V 0 , phase transition from the 3D XY-spin ordered state to the 3D RFIO takes place as σ x is increased. Both the states have the own LROs. On the other hand for larger V 0 , the system has a quasi-LRO for smaller σ x , and the state turns to that of the RFIO with the genuine Z 2 LRO as σ x is increased. In the limit 1 V 0 → 0, the system can be regarded as a classical 2D system. Therefore the study of this limit reproduces the result of the previous study on the classical XY model in 2D with the random external field [3]. There the RFIO forms simply as properties of the lowest energy state.

Robustness of a finite RFIO state
The phase diagram obtained in the previous section, Fig.6, shows that the state of the quasi-LRO changes to the state with the genuine Ising-type RFIO as σ x increases for a large V 0 . This result indicates the robustness of the RFIO state. In order to verify this fact, we study how finite-size systems of the quasi-LRO and also the RFIO change under the updates of the MC simulations. Result of ( S x , S y ) for fixed τ is shown in   Data are plotted for every 10 4 sweeps. Similar result to that of the RFIO in Fig.7 is also obtained in the 2D RFIO state in Fig.6.
From Fig.7, it is obvious that, in the ordinary system without the random external field, the average of the pseudo-spin (magnetization) is nonvanishing but unstable under the MC update, i.e., the orientation of the magnetization fluctuates strongly because of the finiteness (smallness) of the system. On the other hand in the random case, the average is quite stable and stays S y i = ±1, S x i = 0. Then one may wonder that the spin behaves as a classical spin and lost its quantum properties. In order to study it, we measured behavior of ( S x , S y ) as a function of τ . See Fig.8. The results in Fig.8 indicate the following fact. For a small V 0 (i.e., large c τ ), the phase degrees of freedom of the boson operators a i and b i have small quantum fluctuations as their the boson densities fluctuate rather largely. Then, the spin behaves as a classical spin. On the other hand for a large V 0 (small c τ ), the spin behaves as a quantum spin and a superposition of the ↑ spin and ↓ spin is possible. This result indicates that the 2D RFIO state can be used as a quantum qubit in the quantum information device. In Figure 9. (Color online) Methods of making a superposed state of single qubit and an entangled state of two qubits. In the vanishing random-external field, spin is unstable as a result of the quantum fluctuations. By applying a random field suddenly to that state, a superposed state is expected to form. Fig.9, we show a method to make a superposed state | ↑ + | ↓ and also an entangled state of two qubits. Study on quantum superpositions of macroscopically distinct states has the long history [15]. The above study indicates the possibility that a mesoscopic RFIO state is a candidate for the quantum mesoscopic superposed state. In this subsection, we shall study finite-T phase diagram of the random system. In particular, we are interested in how the states with the RFIO evolve as T is increased.

Finite temperature phase diagram
This study is closely related with the stability of the RFIO states investigated in the previous subsection..
In the present MC simulation, the temperature is given as k B T = 1/(N τ ∆τ ). Then, system at low T is realized for sufficiently large N τ . Temperature of the system is increased by decreasing ∆τ for fixed N τ , and therefore the parameters in the action vary as indicated in Eqs. (15). It is obvious that the original 3D system tends to be quasi-1D as ∆τ → small, and then the LROs disappear, which is nothing but a finite-T phase transition.
It is interesting how the ordered states evolve as T is increased. In order to identify the finite-T phase diagram, we measured the internal energy E and the specific heat C as the investigation of the quantum phase transition in the previous section. We also calculated the various correlation functions to identify each phase transition. We show the obtained results and the phase diagram in Fig.10. In Fig.10, for the states of the 2SF for moderate hopping amplitude C a 3 = C b 3 , the specific heat C exhibits two sharp peaks that indicate a second-order phase transition. The correlation functions show that the 2SF state first loses the properties of the BECs and then the pseudo-spin LRO as T is increased. On the other hand for deep 2SF state, a first-oder phase transition takes place from the 2SF to the disordered state directly. This disordered state should be distinguished from the PM state in the low T phase diagram. The latter appears as a result of the competition between V 0 -term and the hopping, whereas the present one comes from the effect of the thermal fluctuations.

Topological excitations in RFIO state
It is interesting to compare topological excitations in the genuine FM state and RFIO state. There are two topological objects that play an important role near phase boundary; (i) vortices of the a and b-atoms and their bound state (ii) domain wall of the relative phase of the BECs of the a and b-atoms It is well known that the a and b-vortices proliferate in the PM state, whereas, in the FM state, the spatial overlap of these two kind of vortex increases as a result of the coherent condensation of the spin operator S x i and/or S y i . Furthermore in a constant external magnetic field h, the Zeeman coupling h· S i generates a linear-potential between the a and b-vortices and "confinement of vortices" takes place. For example h = (h, 0), the Zeeman coupling is given as hS x i = h cos(θ ai − θ bi ), where θ ai (θ bi ) is the phase of the a(b)-atom, and then for a configuration of vortex pair, an extra energy is generated proportional to h·(distance between two vortices in a pair). In two-gap superconductors, a mixing of the two Cooper pairs gives a similar effect to the Zeeman coupling in the FM, and therefore it generates a confinement of vortex [16].
In the system in a random external magnetic field, the Zeeman coupling such asJ x i cos(θ ai − θ bi ) shows up. As seen in the previous section, configurations like (θ ai − θ bi ) ∼ ± π 2 dominates because of the random Zeeman coupling. Therefore it is expected that the interaction between a pair of a-vortex and b-vortex in the RFIO is quantitatively different from that in a constant magnetic field, i.e., a constant Rabi oscillation that prefers configurations with (θ ai − θ bi ) ∼ 0. We investigate this problem in this section. Expected configuration of a vortex pair in the RFIO state is shown in Fig.11. We calculate the local density of vortices V a r and V b r , which is defined as with the vorticity v Ar at site r of the 3D space-time lattice, where A = a, b. Here we have introduce a cutoff and set V A r = 0 if |v Ar | is smaller than 1/2. This cutoff is useful for clarify locations of vortices. From the local vortex density V A r in Eq.(32), we measure the overlap of the vortex configurations of the aand b-atoms by calculating dV in each time slice, which is defined as where N v is the total number of vortex with the Heaviside θ-function, θ(x). We show the calculation of dV and also the total number of vortex and anti-vortex, (N A+ v , N A− v ) in Fig.12. Larger dV means a smaller overlap of the a-and b-vortices.
Transition from the PM state to the RFIO state takes place at C 1 ≃ 0.62. Density of vortex and anti-vortex does not change substantially in the PM and RFIO phases, but the vortex overlap dV changes drastically at the phase boundary. See Fig.12. We also show similar quantities for the nonrandom case with constantJ x i = 0.3. The result indicates that energy of the vortex pair in the RF system is smaller than that in the system of constantJ x i = 0.3. From the above result, it is also expected that a shape of the brick wall between a-and b-vortex pair in the RFIO state is different from that in the constant field, although it is not so easily to observe it by snapshots of the MC simulations. See Fig.11.
Let us turn to the domain wall of the relative phase (θ ar − θ br ). This domain wall is closely related to the "string" connecting a-atom and b-atom vortices. See Fig.11. As we explained above, the calculation of dV indicates that energy of the domain wall is getting smaller as the randomness ofJ x is getting larger from σ x = 0. Then in the case of a large σ x , the pseudo-spin loses its order as a result of a large spatial fluctuation of J x .
To verify that the above expectation is correct, we investigate configurations generated by the boundary condition such that the spins S r on the left spatial boundary have (θ ar − θ br ) = −π/2, whereas on the right spatial boundary (θ ar − θ br ) = π/2. In Fig.13, we show the expectation value of spin S r for various σ x . For the case C 1 = 3.0 and C a 3 = C b 3 = 2.0, the result obviously indicates that, from σ x = 0 to σ x = 1.1, the stiffness of spins is getting stronger as a result of larger fluctuation of the random variableJ x . On the other hand in the case of σ x = 3.0, the pseudo-spins fluctuate rather strongly. Similar behavior is observed in the other cases in Fig.13. This indicates that the strongly fluctuating random-external field destroys the spin order. Then, it is an interesting problem to determine a critical randomness σ c for the order-disorder phase transition observed in the present numerical simulations. Figure 13. (Color online) (Top panel) Expectation value of spin S r , S r , for various σ x with the boundary condition such as ↑ (↓) on the right (left) boundary. Measured length of S r at each site is almost unity. In the region of kinks, length of S r is modified to show the direction of S r clearly. For larger σ x , the random variableJ x fluctuates more strongly. Result obviously shows that in the region of the moderate fluctuation ofJ x , σ x = 0.3 ∼ 1.1, the domain wall is rather thin compared to the case ofJ x = 0 (σ x = 0). However, in the case of σ x = 3.0, the spins fluctuate rather strongly as a result of a large fluctuation ofJ x . System size is L x = 50, L y = 10 and L t = 10. C 1 = 3.0 and C 3 = 2.0. (Middle and bottom panels) Similar behavior is observed near phase boundary for (C 1 , C 3 ) = (0.8, 0.5) and (C 1 , C 3 ) = (1.5, 0.1).

Conclusion
In this paper, we studied effect of a "random external field" on the phase diagram of the bosonic t-J model, properties of the states and the low-energy excitations in the RFIO state. This external field is realized by a random Rabi oscillation between two internal states in an atom induced by a random Raman laser. In the phase diagram of the bosonic t-J model without the random field, there exist ordered states such as the pseudo-spin FM state and the 2SF. We first investigated how the phase diagram is changed as a result of the random field and found that the ordered states move to the states with the RFIO. In the RFIO states, the original U(1) symmetry reduces the Z 2 -Ising type, and therefore low-energy excitations in the RFIO states have different propertied from those in the original ordered states of the t-J model.
By the replica-MFT, we first studied the low-energy properties of the quantum spin system in a random external field, and found that, for a sufficiently strong randomness, there appear the preferred directions of the spin order, which are perpendicular to the applied field. Then, by using the MC simulations, we studied the phase diagram of the effective field theory of the t-J model in applied random external fields, (J x ,J y ). We found that the direction of the spin order is determined by which component of the applied field, (J x ,J y ), is larger. We also studied the finite-T phase diagram and found that the RFIO of the spin survives at intermediate temperatures although the SF is destroyed by the thermal fluctuations.
Finally, physical properties of topological excitations such as the vortex and domain wall were studied. Binding energy of the vortex pair of the a-and b-bosons is smaller in the RFIO compared to that in the genuine t-J model. This means that average distance between a-and b-vortices in a single vortex pair is getting longer as σ x increases. Similarly, the width of the domain wall is thiner in the RFIO state. We hope that the above findings are observed by experiments on cold atomic gases. In near future, we shall report studies on the behavior of vortex lattices that form as a result of the coupling to an artificial external vector potential.