{'fraction_non_alphanumeric': 0.0439722841360597, 'fraction_numerical': 0.010853764899699584, 'mean_word_length': 4.140348692403487, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 3, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Diffusion Monte Carlo methods for Spin-Orbit-Coupled ultracold Bose gases

We present two Diffusion Monte Carlo (DMC) algorithms for systems of ultracold quantum gases featuring synthetic spin-orbit interactions. The first one is a discrete spin generalization of the T- moves spin-orbit DMC, which provides an upper bound to the fixed-phase energy. The second is a spin-integrated DMC method which recovers the fixed-phase property by avoiding the definition of the effective Hamiltonian involved in the T-moves approach. The latter is a more accurate method but it is restricted to spin-independent two-body interactions. We report a comparison between both algorithms for different systems. As a check of the efficiency of both methods, we compare the DMC energies with results obtained with other numerical methods, finding agreement between both estimation


I. INTRODUCTION
The interplay between the electron spin and its momentum, known as spin-orbit coupling (SOC), is an effect of major relevance when studying a wide variety of systems in the field of solid-state physics, such as Majorana fermions [3], spintronic devices [4] or topological insulators [5]. The realization in the last few years of a synthetic SOC interaction in ultracold atomic gases, by exploiting the space-dependent coupling of the atoms with a properly designed configuration of laser beams [6], [7], [8], [9], represents an important achievement. More interestingly, these new realizations allow for a better understanding of the effects induced by the SOC interaction, since ultracold quantum gases are highly controllable and tunable [10]. Ultracold SOC quantum gases have been studied in the dilute regime [11], showing the rise of new exotic phases, such as a spin-polarized plane wave phase and a stripe phase. This stripe phase has been recently observed by Li et al. [12] showing specific properties of a supersolid phase.
Up to now, the theoretical approaches used in the study of SOC gases rely on the mean-field approximation. This theory is expected to be valid when the gas parameter is very small, na 3 ≤ 10 −5 , but beyond this limit one is faced with beyond-mean field terms. A way of surpassing the range of applicability of the mean-field approximation is the use of quantum Monte Carlo (QMC) methods, which are not based on any perturbative scheme. In the present work, we use QMC to study these ultracold atomic gases featuring a synthetic SOC interaction. In particular, we work with the Diffusion Monte Carlo (DMC) method, which is a stochastic method intended for solving the imaginary-time many-body Schrödinger equation. The action of the imaginary-time propaga- * juan.sanchez.baena@upc.edu † jordi.boronat@upc.edu ‡ ferran.mazzanti@upc.edu tor exp −τĤ is implemented as a set of transformations to a list of points in coordinate space (commonly called walkers) that represent statistically the wave function. In the limit τ → ∞, the ground state dominates while excited-state contributions are exponentially damped, providing an exact estimate of the ground-state energy and of any observable commuting withĤ. If the ground state of the system of interest is complex (which is the case when the SOC term is present), it is necessary to invoke the fixed-phase approximation (FPA), which provides an upper bound to the ground-state energy. Previous DMC calculations with SOC terms in the Hamiltonian have been carried out in the study of electronic structures [1,13], quantum dots in semiconductors [14], and repulsive Fermi gases [15]. A DMC method incorporating the SOC terms that arise in electronic systems has already been developed [1]. In this method, the authors implement the spin-orbit term of the propagator through the use of the T-moves technique [2]. They also use a regularized, continuous representation of the spin degrees of freedom. In order to control the sign problem that the SOC terms introduce in the propagator, the authors of Ref. [1] define an effective Hamiltonian in such a way that the propagator becomes positive-definite. It can be shown that the estimations obtained with this effective Hamiltonian yield an upper bound to the fixedphase energy [16]. In the present paper, we adapt the T-moves DMC algorithm of Ref. [1] to the usual, discrete representation of the spin, and show how to treat the synthetic SOC present in ultracold quantum gases. We also introduce a different method for treating the SOC terms of the propagator, loosely based on Ref. [17], which consists on propagating the wave function integrated over all spin configurations. In doing so, we avoid almost completely the sign problem induced by SOC terms, meaning that no effective Hamiltonian needs to be defined. This paper is organized as follows. In Sec. II, we discuss the form of the Hamiltonian as well as several kinds of spin-orbit couplings of interest in the field of cold Bose gases. The reduced units used in this work are introduced in Sec. II A. In Sec. III we present the details concerning the Spin-integrated DMC. We derive the fixed-phase, spin-orbit propagator to first order and elaborate on its implementation within the DMC framework. We also discuss how to estimate the energy and provide a scheme of the DMC algorithm. In Sec. IV, we show how to implement discrete spin sampling within the T-moves DMC , as well as how to implement the SOC terms introduced in Sec. II. In Sec. V, we compare both DMC methods in one and two-body problems (V A) and in some manybody cases (V B). Finally, in Sec. VI, we summarize the main conclusions of our work.

II. HAMILTONIAN
The system studied in this work is formed by an ultracold gas of N bosons of mass M with pseudo-spin 1/2 under the effect of synthetic spin-orbit coupling [10]. The generic form of the Hamiltonian is: withV 1b k andV 2b momentum independent, local, oneand two-body interactions, respectively. Notice thatV 2b can depend on the spin configuration. In much the same way,Ŵ SOC k stands for a one-body, momentum and spindependent potential. The ones considered in this work are the Rashba, Weyl and Raman interactions given bŷ withP α k the α-component of the momentum operator of particle k,σ x,y,z k the Pauli matrices associated to particle k, Ω the Rabi frequency, and λ α (α = { Rs, We, Rm}) the strength of the corresponding SOC interaction. The general form of the two-body potential is: where s k , s l assign values ±1 to the z-component of the spin of particles k and l, while V 2b s k ,s l (r kl ) is a central, short-ranged potential that can be different for the different channels corresponding to s k and s l . In the numerical examples of Sec. V we use a soft-core force, defined by: If the two-body interaction is taken to be spinindependent, V 0 (s k , s l ) = V 0 and R 0 (s k , s l ) = R 0 .
The one-body potential used in some of the calculations below is:

A. Reduced units for the different kinds of SOC interactions
Due to the different spin dependence, we use different length and energy scales in each case. These are the following: for the Rashba interaction, we set the length and energy units to while for the Raman interaction with η Rm a dimensionless scaling factor that we vary depending on the density. Finally, for the Weyl Hamiltonian we use In terms of these, the interactions read where all quantities are dimensionless. The same applies to the soft-core potential and harmonic trap of Eqs. (5) and (6).

III. THE SPIN-INTEGRATED DMC (SIDMC) METHOD
A. The spin-orbit propagator in the fixed-phase approximation In this section we derive a suitable form of the propagator required to simulate spin-orbit problems, under the assumption that the two-body interactionV 2b is spin-independent. The imaginary time evolution of state |Ψ(τ ) is given by Projecting on R , S | and introducing an identity, Eq. 13 can be written as: where R and S stand for the position and spin coordinates of the N particles. For the sake of clarity, we also defineĤ Up to O(∆τ ), Eq. (14) can be written as: where the term corresponding toĤ 0 in the splitting of H in the propagator is spin-independent. In this way, the propagator reads: This propagator can have complex contributions coming from the Pauli matrices appearing in the spin-orbit interaction, and therefore sampling it is not possible. In order to bypass this problem, we resort to the fixed-phase approximation [1] where all quantities involved are real.
Knowing the general expression of the propagator written above, we can deduce its reduction to the fixedphase approximation. This can be done comparing the imaginary-time Schrödinger equation for the wave function and for its magnitude, which is the main quantity of interest in the FPA. For the full wave function, one has while for the magnitude ρ( R, S) of ψ( R, S) the equation reads where ψ( R, S, τ ) = ρ( R, S, τ ) exp iΦ( R, S, τ ) and R, S|ŵ Re | R , S = Re R, S|Ŵ spin | R , S e iΦ( R , S ,τ ) e iΦ( R, S,τ ) .
(22) In the FPA, Φ( R, S, τ ) is independent of τ andV Φ = N k=1 |∇ k Φ( R, S, τ )| 2 becomes a local interaction in positions and spins. Equations (19) and (20) have a similar structure, and thus comparing terms in each, we can get the FPA form of the propagator in Eq. (18): withĤ the fixed-phase Hamiltonian. In the FPA one has to impose a certain form for the phase. In this work we impose it to be the sum of one-body terms Due to the form of the spin-orbit potential, we can evaluate the integral in Eq. 23. For the Raman SOC of Eq. 2, the matrix element of the spin-dependent part of the potential is where Since the spinless part of the propagator is given by [18], one has with V 0 the spinless part of the potential entering inĤ 0 and E s the common energy shift used in the DMC algorithm. Up to O(∆τ ), the integral in Eq. (23) yields For the Rashba and Weyl SOC interactions, a similar procedure has to be carried out. However, one has to expand the element R , S | exp −∆τ ŵ Re +V Φ | R , S in Eq. (23) up to order ∆τ 2 . This is because the terms originated from the matrix element ofŵ Re are proportional to ξ k = r k − r k , and thus, the elements arising fromŵ 2 Re generate contributions of order ξ 2 k and ξ k ξ l . Since ξ k represents the displacement of particle k due to the standard DMC Gauss-Drift-Branching (GDB) process, this quantity is of O √ ∆τ . However, in the numerical experiments conducted, we have not found a significant impact on the results when these terms are dropped.
Following with the derivation of the propagator in Eq. (29), we define a new operatorÔ as while, up to O(∆τ ), Eq. (29) can be rewritten as: where we have used the approximation (1 − ∆t x i ) ≈ (1 − ∆tx i ) which is exact to order ∆t. In this way, the matrix element of the new operatorÔ becomes the product of matrix elements of single-particle operatorŝ O k , as shown in the expression above.
Note that, for the Rashba and Weyl SOCs, the matrix elements s k |Ô k |s k depend both on r k and r k . For the sake of simplicity, in the following we omit the r k and r k labels. The imaginary time evolution equation for the magnitude of the wave function, within the fixed-phase approximation and to order ∆t, becomes However, in DMC simulations the object that is propagated is f ( R, S, τ ) = ρ( R, S, τ )ρ T ( R, S), with ρ T ( R, S) the magnitude of a given importance sampling trial function. From Eq. (32) one readily sees that In order to implement this equation, we need the propagator to be positive-definite. However, due to the spinorbit coupling, the matrix elements of the propagator do not fulfill this condition. Despite this, if we propagate the spin-integrated form of the magnitude of the importance sampling function f of Eq. (33), this problem is greatly reduced. Therefore, we propagate the quantity In order to progress, we impose the magnitude of the trial wave function to be spin-independent i.e. ρ T ( R, S) = ρ T ( R). After j time steps, one gets where R (n) are the position coordinates of the walker, and s (n) k the spin of particle k of that walker, both at iteration n. We can understand this expression in a simple way. The last two pieces correspond to a standard GDB DMC process [18] for the spinless part of the Hamiltonian. On the other hand, the first part, incorporating the spin-dependent terms, can be implemented through a secondary branching process. This one must fulfill that, after j iterations, the weight carried by a given walker is given by corresponding to the first term in Eq. (35). This is fulfilled by performing the secondary branching at iteration j using the weight with the initial condition w(0) = 1. It can be shown that w(j) can be easily computed as in terms of the spin weight factors where |↑ and |↓ stand for |s = 1 and |s = −1 , respectively. In this way, in the proposed method each walker carries the evolution of both c + and c − for every particle, instead of explicit spin variables.
Notice that the weights in Eq. (38) are the products of one-body terms, due to the one-body nature of the spindependent part of the interaction. Therefore, for each k, we can write Equation (36) can be expressed in terms of η + k (j) and η − k (j), and we can rewrite Eq. (35) as It can be shown that the marginal spin integrated magnitude of the importance sampling function of Eq. (33) reads where S N −1 in the second term means summing over the spins of all particles but the k-th one. This quantity is relevant in the energy estimation, described in the next Section.

B. Energy estimation
We show in this Section how to estimate the energy of a many-body system under SOC interactions using the method introduced in the previous Sections, although it can be easily extended to estimate any other quantity. The DMC energy estimator in the FPA at iteration j is given by: withĤ FP defined in Eq. (24). The local energy is, therefore, which, as it can be seen, depends on R (j) and S, so that . We can split it in two parts corresponding the the spin-independent and spindependent contributions, respectively. The spinindependent part can be expressed in the form while with L,S,l the one-body contribution to the spindependent local energy corresponding to particle l (recall thatŵ Re +V Φ is a one-body operator). With all these definitions, Eq. (43) becomes The term E DMC,0 (j) contains all the spin-independent contributions, and can be written as with F ( R, τ ) defined in Eq. (34). This part of the energy is evaluated as usual in DMC, i.e.
where N w is the total number of walkers in the simulation, and i w specifies the walker index. In much the same way withF ( R (j) , s l , j∆τ ) defined in Eq. (42). Therefore, we need to be able to sampleF ( R (j) , s l , j∆τ ) in order to evaluate E DMC,S (j). This can be done by estimating E DMC,S (j) as This expression ensures that each local energy contribution is averaged with an effective weight given by which is the one associated toF ( R (j) , s l = ±1, j∆τ ) in Eq. (42). It is important to realize that Eq. (52) can be used to estimate the expectation value of any quantity that depends on the spin through one-body terms only, replacing the L,S ( R, ±1) terms with the corresponding operators.

C. The SIDMC algorithm
In this section we present a scheme of the Spinintegrated DMC algorithm. In the present method, a walker is represented by the set of quantities Particle positions are initialized as usual in Monte Carlo simulations, while spin weight factors c ± k must be initialized to one in the first iteration The first step in each iteration of the algorithm is to perform a standard GDB process using the spinless part of the HamiltonianĤ 0 and ρ T ( R). Next, one has to update the c ± k coefficients according to the expression which yields the new coefficients at iteration j + 1 from the known ones at iteration j. Notice that, in this expression, O k is the 2 × 2 matrix of Eq. (39). Once with these coefficients, one can obtain w(j + 1) according to and from here, the secondary branching factor, Notice this weight is different for each walker, so in fact B = B iw with i w the walker index.
In practice, it may happen that, along the simulation, the absolute value of the c ± k (j) coefficients keeps increasing unboundedly. However, the ratio of w's in this equation is always finite. On the other hand, it is better to use a mixed-branching strategy with the B(j + 1) terms, where walkers acquire a weight that is being updated along each block of iterations. The accumulated weight B iw at the end of the block is equal to the product of the weights at each iteration, for each walker. Once the block is finished, these weights are used to replicate the list of walkers.
In DMC simulations, the weight of the walkers is divided by a constant (equal to e E T ∆τ with E T the threshold energy and ∆τ the time step) when performing the replication process [18]. One has to perform an equivalent renormalization with the secondary branching, while in this case the normalization constant can be computed in two ways. One way is to use the average over the final number of walkers of the accumulated B of the previous block. Another way is to use the B coefficients of the current block, accumulated over the previous iterations and averaged over the number of walkers. The best strategy is determined by the SOC model at hand, with the first choice being more suitable for the Raman interaction, and the latter performing better with the Weyl and Rashba models.
The energy at iteration i inside a block is estimated as with E (iw) L,0 ( R (j) ) and ε (iw) L,S ( R (j) ) given in Eqs. (46) and (52). In this expression the sum is over the complete set of N w walkers, obtained after the standard GDB process associated to the spinless part of the Hamiltonian.
In this way, the expression implicitly includes the weighting of the standard branching. This equation represents the generalization of Eqs. (50) and (52) for the mixedbranching case.
An important remark concerning the secondary branching is that B(j + 1) in Eq. (58) is not positive definite. However, the fraction of walkers which generate a change in sign is tiny, and thus walkers that produce this effect can be safely discarded. To quantify that, we monitor the quantity with N e and N b the number of eliminated walkers and the number of iterations per block, and N w the average number of walkers of the block. Our numerical results show that χ depends slightly on the value of the parameters chosen for the simulation, but it is always of the order of 10 −3 or smaller.

IV. DISCRETE SPIN T-MOVES DMC (DTDMC)
In this section we adapt the continuous spin T-moves method of Ref.
[1] to a system of discrete spins under the SOC interactions analyzed in this work. In the following, we assume the two-body interaction is spin-dependent, with (possibly) different contributions in each channel. In this method the walkers carry explicit spin variables together with the particle positions.

A. Formalism
In order to derive the alternative algorithm, one has to go back to the beginning and work out the propagator in Eq. (23), which we split in a different way rearranging terms as follows We can introduce the importance sampling function inside this expression and write To order O(∆τ ), the first term inside the integral becomes However, for any kind of spin-orbit coupling the matrix element R , S |ŵ Re | R , S is not always negative, and thus Eq. (65) can not be interpreted as a probability distribution. In order to bypass this limitation and in the spirit of Refs. [1,2,16], we define an effective Hamiltonian that replaces the original one, and that leads to a variational upper bound to the fixed phase energy of the original Hamiltonian. We thus writê where the sumŵ eff Re,A +ŵ eff Re,B is an approximation to the originalŵ Re of Eq. (22). This approximation is built such that the local energy ofĤ FP eff andĤ FP are equal when they act on the magnitude of the trail wave function. The matrix elements of these terms are given by with the transition coefficients where in the last expression, the summation and the integration are restricted to those values that satisfy the condition T > 0. Using these definitions we avoid non-local matrix elements producing negative transition probabilities. Notice also that the effective Hamiltonian depends on the magnitude of the trial wave function, which means that the the energy obtained depends on its choice. We showcase this effect in Sec. V. The fixed-phase propagator for the effective Hamiltonian, with importance sampling, is thus: Since this propagator is positive-definite, we can now interpret it as a probability distribution. Therefore, one can sample from it. This can be implemented performing initially a GDB of the exp −∆τ (Ĥ 1 +ŵ eff Re,B ) part, with a branching factor that, according to Ref. [2], reads which generates the displacement R → R . In a second step, one performs a transition ( R , S) → ( R , S ) given by the probability where Despite the sum in Eq. (73) involves the 2 N spin configurations, which sounds prohibitive for large N , it must be kept in mind that only one-body operators are involved and the expression is greatly simplified.

B. Application to synthetic SOC in ultracold gases
In this section we particularize the results of the previous formalism to the Weyl SOC interaction. The procedure is analogous for the Rashba and Raman potentials. We start evaluating the matrix elements ofŵ Re , which are given by with φ k the single-particle phase of Eq.( 25). In this expression we have omitted the last term of Eq. (3) as it is a constant contribution that represents a shift of the total energy only. In order to construct the effective Hamiltonian, we must evaluate the matrix elements of w Re to check their sign. However, given any set of coordinates r k , r k , terms of the form d dx k (δ(x k − x k )) are in general problematic. In order to preserve the upper bound property of the effective Hamiltonian, we adopt the (apparently rude) prescription (76) with a small parameter. This is equivalent to replacing the momentum operator witĥ while both expressions coincide to order . Notice that, in this form, the resulting operator is still hermitian, and for → 0, the energy is preserved. With this substitution, w Re is replaced by a new operatorŵ Re, , whose matrix elements are the same as in Eq.
We discuss in this Section a scheme of the DTDMC algorithm to better understand its practical implementation. A walker at iteration j is described by with s k = ±1 the z-component of the spin of particle k and subindexes and superindexes standing particles and iterations, respectively. The initial condition for the position and spin coordinates is generally obtained through the sampling of the trial wave function using the Metropolis algorithm.
The first step to be implemented at each iteration is a GDB process with the branching factor given by Eq. (71), which produces a spatial translation R (j) → R (j) A . After this, we need to sample the part of the propagator which depends on the effective potentialŵ Re, ,A . In this second step, a transition ( R

(j)
A , S = S (j) , R = R (j+1) and S = S (j+1) . As an example, we explicitly report how this evolution is carried out for the Weyl SOC case. A possible transition probability is: Notice that the terms appearing in P ( R, S → R , S ) are different for each walker and each iteration. In general, one has to keep here only those terms of Eq. (75) (after the substitution of Eqs. (76) and (77)) that are strictly negative. This total transition probability is the sum of different transition probabilities P (m) t,k , so it can be written as The probabilities P (m) t,k depend on the coordinates of all particles but we only make explicit the dependence on the coordinates that change under each transition for the sake of simplicity. Notice that in this example there are 4N + 1 possible transitions. We define the cumulative distribution vector as Notice that v c (i c ) ∈ (0, 1] ∀i c . To sample this discrete probability distribution function we follow the standard procedure: we generate a random number ξ ∈ [0, 1] and select the component Finally, we perform the transition associated to the t,k , the spin of particle k flips and its coordinates are modified according to x k = x k , y k = y k + , z k = z k , while the rest of the system is left unchanged.

V. RESULTS
We report in this Section results for the energy in different systems for both the SIDMC and DTDMC methods. In Sec. V A, we show the energy of a few one-body and two-body problems, while in Sec. V B, we report results for the energy of a few many-body systems, both in the mean-field regime and out of it. As a check of validity of the two DMC algorithms for SOC systems, we compare the DMC estimations with energies obtained from the imaginary-time evolution of the Schrödinger equation (one and two-body cases) and the Gross-Pitaevskii equation (many-body in the dilute regime). We also comment on the technical issues mentioned in Secs. III C and IV C, mainly the elimination of walkers in SIDMC and the influence of the parameter in DTDMC, as well as the dependence of the energy estimation on the time step. In all cases, the parameters of the Hamiltonian and the trial wave function are reported in reduced units (see Sec. II A).

A. One and two-body problems
In this Section, we report DMC results for the energy corresponding to four different physical situations: a three-dimensional (3D) one-body system with Weyl SOC, a 3D one-body system with Raman SOC, and two interacting two-dimensional (2D) two-body systems with Rashba SOC, one featuring a spin-independent two-body interaction and another with a spin-dependent one. All systems are harmonically confined. We summarize our results in Table I, which includes the DMC energies obtained with both algorithms together with the imaginary time evolution (ITE) estimates, both for the fixedphase Hamiltonian (Eq. (24)) and the fixed-phase, effective Hamiltonian (Eq. (66)). All SIDMC energies are obtained by performing several simulations, changing the parameter ∆τ , and then extrapolating the energy to the limit ∆τ → 0. In the Weyl and Rashba cases with DT-DMC, one must carry out several calculations changing ∆τ and and then extrapolate to the limits ∆τ → 0, → 0, and ∆τ → 0. We discuss below how to perform the triple limit involving ∆τ , , and ∆τ . This setup is not necessary in the Raman calculations since the SOC part of the propagator scales as O(N ∆τ ) if is sufficiently small.
The trial wave function for each Hamiltonian is important because it fixes the phase and, in all cases, reduces the variance via importance sampling. In the problem of Raman SOC and DTDMC the trial wave function that we have used is with µ = 1 2 acos k ηRm , k the reduced momentum and ω the harmonic oscillator strength. In these expressions, {k, C 1 , C 2 } are taken as variational parameters. The SOC term of the trial wave function is of the same form as the one used in Ref. [10]. Since the magnitude of the trial wave function must be independent of the spin in SIDMC, we have used with B c another variational parameter.
Concerning the Weyl model, the adopted trial wave function for DTDMC is where θ k and φ k are the polar and azimuthal angles of the momentum vector k, respectively. The adopted magnitude of the trial wave function for the SIDMC case is Finally, the trial wave function used in the DTDMC two-body Rashba simulations is with φ k the angle of the momentum vector in polar coordinates. In this expression, ρ T,2b ( r 1 , r 2 ) is the exact solution of the two-body interacting problem at low momentum (k 2b ∼ 10 −2 ) (without SOC) corresponding to the soft-sphere potential of Eq. (5), with parameters This choice makes the two-body trial wave function spinindependent for simplicity. We use the same choice for the SIDMC simulations. The time step is ∆τ ∼ O 10 −3 in DTDMC simulations while it is ∆τ ∼ O 10 −2 in the SIDMC ones. The average number of walkers is kept stable along the simulations, and it is fixed to a value between 2000 and 3000, depending on the case. The parameter of DT-DMC is fixed as = 100∆τ in the Raman calculation and as = 200∆τ in the Rashba and Weyl cases. In the Weyl SIDMC calculations, the secondary branching weights w(j) are accumulated along blocks of N b = 10 iterations. The ratio of eliminated walkers is χ < 0.001. In the Rashba cases, we have N b = 50 and χ < 0.002. Finally, for the Raman problem we have N b = 10 and χ = 0 (see Sec. III C).
In Fig. 1, we show the energy as a function of the imaginary-time step for the two-body Rashba calculations. We can clearly see a linear dependence of the energy with the time step, as it corresponds to a linear approximation to the exact propagator. In the DTDMC method, as stated previously, three limits have to be satisfied in order to obtain the estimation of the energy: ∆τ → 0, → 0, and ∆τ → 0. The extrapolations according to these limits can be performed in several ways. Here, we present two of them. Method 1 consists on performing N sets sets of N sim simulations making ∆τ → 0, → 0, with ∆τ 1 fixed. After this, one ends up with N sets estimations of the energy, each one associated to a given ∆τ value. Finally, one retains the estimation associated to the lowest ∆τ value. Method 2 consists on performing N sets sets of N sim simulations making ∆τ → 0, ∆τ → 0, with fixed. After this, one ends up with N sets estimations of the energy, each one associated to a given value. Finally, one then takes the extrapolation of these estimations in the limit → 0.
In Figs. 2 and 3, we show the estimations obtained using Method 1 and Method 2, respectively, for the onebody system with Weyl SOC and a harmonic trap. As we can see, the dependence of the energy extrapolations with respect to ∆τ is much weaker than their dependence on . Therefore, Method 1 is preferred and is the one that we have used to provide the T-moves energy. We can also see from the figure that the dependence of the energy with respect to ∆τ , when or ∆τ are fixed, is linear in both cases. This is because the non-SOC terms of the propagator are exact up to O(∆τ ) while the SOC terms are exact up to O ∆τ . For all the chosen values of , the conditions in Eq. (78) are satisfied, with the r.h.s being 10 −2 . Also, the condition in Eq. (79) is satisfied since the difference between the r.h.s. and the l.h.s. is at most a 3% of the SOC local energy contribution.
From Tab. I, we can see that both DMC methods provide energies that agree with the result of the imaginarytime evolution within a 2% error. We can also see that SIDMC provides lower energies than DTDMC. This is due to the fixed-phase nature of the energies obtained with SIDMC, which does not require to use an effective Hamiltonian as DTDMC. We can see that this effect is enhanced in the harmonically trapped systems featuring Rashba and Weyl SOCs. For the cases with two-body spin-dependent interactions, only T-moves results are re- ported, since SIDMC can not deal with these kind of potentials. It must be remarked that, while in the T-moves calculations we perform the triple extrapolation ∆τ → 0, → 0, and ∆τ → 0, calculations with ITE are performed at a fixed ( = 0.1 and = 0.3 in the Weyl and Rashba cases, respectively). This is due to the computational cost of decreasing when discretizing the Schrödinger equation in the position representation, since is taken as the point-to-point distance of the mesh. In order to check that both DTDMC and ITE give compatible estimates, we also provide in Tab. I DMC energies corresponding to a fixed . This is not necessary in the Raman case since the Raman Hamiltonian is independent of if this parameter is sufficiently small, as mentioned previously. Notice also that the errors corresponding to the ITE results in the two-body 2D Rashba cases are larger than the ones in the 3D one-body Raman and Weyl cases. This is due to the higher number of dimensions that must be discretized in the latter case.

B. Many-body calculations
We report in this Section the DMC energies corresponding to the many-body Raman and Weyl SOC Hamiltonians. We first focus on the dilute regime with a finite number of particles imposing periodic boundary conditions (PBC). We compare the DMC energy estimations with energies obtained by solving the imaginary time Gross-Pitaevskii equation (GPE), both for the fixedphase Hamiltonian (Eq. (24)) and the fixed-phase, effective Hamiltonian (Eq. (66)). In the case of Rashba SOC, we do not know the scattering length of the complete interaction, and thus a direct comparison to GPE is not possible. Finally, we compare the energy estimations of both DMC methods out of the dilute regime. Table II reports the DMC energy per particle together with the corresponding Gross-Pitaevskii energy per particle, for four different physical systems: Raman SOC and Weyl SOC, both with spin-independent and spindependent two-body interactions. Moreover, we include the T-moves energy per particle using two different trial wave functions in the two-body spin-independent Weyl case in order to showcase the variational dependence of this method with respect to the magnitude of the trial wave function.

Dilute regime
For the GPE calculations involving Raman or Weyl SOCs, we use the free-space scattering length, i.e., the scattering length obtained for the Hamiltonian removing the SOC terms [19,20].
In all cases, the trial wave function is of the form with ρ T,2b (r ij ) = ρ T ,2b (rij )+ρ T ,2b (L−rij ) 2ρ T ,2b (L/2) and r ij = | r j − r j |. The function ρ T,2b (r ij ) is the magnitude of a spin-independent two-body trial wave function analogous to the one presented in Sec. V A (here k 2b ∼ 10 −6 ). The magnitude of the one-body terms for the T-moves "Raman 2-b no spin" and "Raman 2b spin" cases are given in Eqs. The average number of walkers is set to N w = 1, 000 and the time step is ∆τ ∼ O 10 −3 . The parameter of DTDMC is fixed as = 100∆τ . All the used values of satisfy the condition of Eq. (79), with a discrepancy between the r.h.s. and the l.h.s. of at most ∼ 1%. Also, the r.h.s of both expressions in Eq. (78) equals 0.08 at most, which implies that the maximum error in the approximation to the propagator is e max ∼ e 0.08 − (1 + 0.08) 0.0033. In the Weyl SIDMC calculations, the length of a simulation block is set to N b = 10. The ratio of eliminated walkers is χ < 0.0002. For the Raman calculations, we have N b = 10 and χ = 0 (see Sec. III C).
We can see from Table II that the DMC energies agree with the GPE calculations up to a ∼ 1%. As in the previous Section, for the spin-dependent two-body cases only T-moves results are reported, since the SIDMC method can not solve two-body spin-dependent interactions. We can also see from the two-body spin-independent cases that DTDMC is able to recover almost completely the fixed-phase energy, although we know it always provides an upper bound to it. On the other hand, SIDMC recovers the complete fixed-phase energy. The DTDMC Weyl two-body spin-independent calculations illustrate the variational property with respect to the magnitude of the trial wave function of this method. Notice that two different magnitudes ("Weyl 2-b no spin" and "Weyl 2-b no spin" cases) provide two different energy estimations.

Beyond the dilute regime
In this Section we compare the performance of the two DMC algorithms discussed in several homogeneous many-body systems, beyond the dilute regime. We analyze a few systems featuring Raman and Weyl SOCs using periodic boundary conditions, and a two-body spinindependent interaction. We show again an example of the variation of the T-moves energy when the magnitude of the trial wave function is changed. We also provide DTDMC energy estimations of systems under Raman and Weyl SOCs with a spin-dependent two-body interaction. Finally, we compare both DMC estimations in a many-body harmonically confined system with Weyl SOC. Results are presented in Table III.
The general form of the trial wave function is given in Eq. (106). The T-moves calculations corresponding to the cases "Raman PBC 2-b no spin", "Raman PBC 2-b spin trial 1" and "Raman PBC 2-b spin trial 2" use the one-body terms of Eqs. (89) and (90), while for the SIDMC "Raman PBC 2-b no spin" calculation Eq. (93) has been used. For DTDMC corresponding to the cases "Weyl PBC 2-b no spin", "Weyl PBC 2-b spin", and "Weyl HO 2-b no spin" we use the expressions in Eqs. (94) and (95) while for the DTDMC "Weyl PBC 2-b no spin trial 2" case we use ρ T ( r, s = +1) = γ (108) This form helps us to illustrate the variational property of the T-moves method with respect to the magnitude of the trial wave function. The SIDMC "Weyl PBC 2-b no spin" and "Weyl HO 2-b no spin" calculations use the expressions in Eq. (98). As in the previous Section, the trial phases for each case are given in Eqs. (91), (92), (96), and (97).
In the two-body spin-independent calculations, the two-body trial terms in all PBC cases are the same as in Sec. V B 1. Concerning the two-body spin-dependent calculations, we report the energy in the Weyl case using a spin-independent two-body correlation factor analogous to the one in Sec. V B 1. In the Raman case, though, we compare the energy estimated using a spin-independent two-body factor with that estimated using a spin-dependent one, again with the same form as in Sec. V B 1. Finally, in the "Weyl HO 2-body no spin" case we set ρ T,2b (r ij ) = ρ T,2b (r ij ) in Eq. (107) because we do not impose PBC.
The average number of walkers is set to N w = 1000, the time step ∆τ ∈ (10 −4 , 10 −3 ), and the DTDMC parameter is fixed such that ∆τ ∈ (100, 400) for Weyl and ∆τ = 10 for Raman. All the used values of satisfy the condition in Eq. (79), with a discrepancy between the r.h.s. and the l.h.s. of at most 3%. Also, the r.h.s of both expressions in Eq. (78) equals 0.3 at most, which implies that the maximum error in the approximation to the propagator is e max ∼ e 0.3 − (1 + 0.3) 0.05. In the Weyl PBC SIDMC calculations the length of a simulation block is set to N b = 10. The ratio of eliminated walkers is χ < 0.006. The harmonically trapped Weyl simulations share the same parameters except for the ratio of eliminated walkers, χ < 0.001. For the Raman calculations one has N b = 10 and χ = 0 (see Sec. III C).
In Fig. 4, we show the energy dependence on the imaginary time-step corresponding to the SIDMC simulations of trapped Weyl gases. We can see in the Figure the linear dependence of the energy with respect to ∆τ . In Figs. 5 and 6, we show DTDMC results for the two methods mentioned in Sec. V A to estimate the triple limit ∆τ → 0, → 0, and ∆τ → 0. The observed behavior is consistent with the previous results obtained in the one-body case.
In Table III, we report the DMC energies for the analyzed cases. From these results, we can see that DTDMC is able to almost exactly recover the fixed-phase energy of the bulk gases. In the trapped Weyl gas, the difference with respect to the fixed-phase energy obtained with SIDMC is larger. We can also see how the improvement of the magnitude of the trial wave function in the twobody spin-independent PBC Weyl simulation produces better energies as a consequence of the variational property of the DTDMC method. Finally, our results show that the spin-dependent two-body trial correlation factor does not make any significant difference in the two-body spin-dependent PBC Raman simulation.

VI. CONCLUSIONS
In this paper, we discuss two different Diffusion Monte Carlo methods (DTDMC and SIDMC) that are able to deal with many-body systems of ultracold quantum gases featuring synthetic Spin-Orbit Coupling. DTDMC is an extended version of the method of Refs. [1] and [2] to the relevant SOC interactions in the field of ultracold gases, but with discrete spins. This method relies on the introduction of an effective Hamiltonian and provides an upper bound to the fixed-phase energy of the system. On the contrary, the SIDMC method is able to avoid this issue by propagating the spin-integrated probability density, providing exact fixed-phase estimations. However, SIDMC is not able to deal with spin-dependent two-body interactions and requires the use of spin-independent trial wave functions.
We have described the formalism of both methods in detail, together with a scheme of both algorithms for future applications. We have reported the energy estimation of several few-body systems, featuring three different kinds of SOC interactions. We have compared these results with energies obtained by propagating the Schrödinger equation in imaginary time, finding good agreement between both estimations. We have also performed simulations of many-body systems in the dilute regime and have recovered the energies obtained by solving the imaginary time Gross-Pitaevskii equation with discrepancies of at most ∼ 1%. Finally, we have compared both algorithms beyond the dilute regime, showing that the DTDMC method is able to recover the fixedphase energy almost completely in the PBC cases. We hope that these methods can be used to explore the physics of SOC systems beyond the mean field, dilute regime.