{'fraction_non_alphanumeric': 0.03418680475642501, 'fraction_numerical': 0.016398158803222096, 'mean_word_length': 4.253652392947103, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Decay modes of two repulsively interacting bosons

We study the decay of two repulsively interacting bosons tunneling through a delta potential barrier by direct numerical solution of the time-dependent Schr\"odinger equation. The solutions are analyzed according to the regions of particle presence: both particles inside the trap (in-in), one particle in and one particle out (in-out), and both particles outside (out-out). It is shown that the in-in probability is dominated by exponential decay, and its decay rate is predicted very well from outgoing boundary conditions. Up to a certain range of interaction strength the decay of in-out probability is dominated by the single particle decay mode. The decay mechanisms are adequately described by simple models.


I. INTRODUCTION
The decay of a particle by tunneling through a potential barrier into a continuum is a fundamental and unique phenomenon in quantum mechanics. The tunneling of multiparticle systems is just as important and presents one of the places where the understanding of macroscopic quantum phenomena can start [1]. The tunneling and decay of Bose-Einstein condensates (BECs) are attractive subjects of study [2], since the BEC is a unique state of matter where quantum mechanical features are manifested at the macroscopic level. After BECs were first realized experimentally in dilute atomic gases [3], a huge amount of related research followed. Ultra-cold atoms are usually trapped in a finite potential well and the decay by tunneling into a continuum is an existing and potentially desirable possibility. In this context it was realized that understanding the decay dynamics by tunneling is a very important task [4,5].
In most cases BECs have thousands to millions of particles and the dynamics is adequately described by the Gross-Pitaevskii (GP) equation [6,7] of mean-field theory [8], a nonlinear Schrödinger equation. The GP equation governs the time evolution of phase and particle number density of an essentially fully Bose-condensed system. With many works on the mean-field description of BEC tunneling [9][10][11][12][13], it is remarkable that there is still a discussion, both about the technical implementation [12] and the correct formulation of mean-field theory related to the decay problem [10]. It is thus desirable to obtain a detailed understanding of the microscopic physics of multi-particle decay.
The cases of stronger interactions or fewer particle numbers are also important, where the GP equation is less accurate. In the few boson regime, correlated decay of particles was observed and studied both experimentally and theoretically [14,15]. The particle correlation in the decayed wave is important in relation to the atom laser [16]. For strongly-interacting bosons in a one-dimensional trap Bose-condensation is not relevant but the gas acquires properties related to fermionic systems [17]. In the Tonks-Girardeau limit of infinite interactions the few boson decay problem was treated analytically [18], and numerical simulation have addressed the crossover for finite interactions from a harmonic trap with up to four bosons [19]. The analytic treatment of few boson decay with finite interaction strength remains a difficult task.
In this paper, we approach this problem by both numerically and analytically. We study the simplest case of two repulsively interacting bosons in a potential trap in one dimension. The time evolution of the decay is obtained from first principles by solving the timedependent Schrödinger equation numerically. Then, it is compared to approximate analytic methods, starting from the exact solutions of local spatial regions. The decay phenomena are investigated for a wide range of interaction strength, from zero to very strong repulsion.
The analytic model predicts exponential decay mode of the interacting system, which is in very good agreement with our numerical simulation. Also, the decay of the total particle number is well explained with a simple theoretical model.

II. THE MODEL HAMILTONIAN
We choose a model Hamiltonian for the two interacting boson decay. Considering the kinetic energy, external potentialṼ ex for trapping and interactionŨ between bosons, the total Hamiltonian with two identical bosons is written as The external potential isṼ and acts as potential trap by a delta barrier at position L. This choice of external potential has some advantages in that the delta barrier has zero width so the consideration of decay process inside the barrier is not needed. Also, the analytical treatment of the decay process is simplified. In single particle case we found that the exponential decay mode dominates and non-exponential features are strongly suppressed compared to a finite-width barrier case.
Computationally, the narrow width of the delta function makes the Hamiltonian matrix more sparse, which makes the problem tractable.
Considering only s-wave scattering [2], the interaction potential between particle 1 and particle 2 is given asŨ whereg is a coupling constant andx 1 andx 2 are the positions of each boson, respectively.
To simplify the analysis and compare the result with external parameters, we introduce dimensionless units. The new length unit x is defined as x ≡x/L. The Hamiltonian is rewritten asH Dividing both sides by 2 /(mL 2 ), we get the rescaled, dimensionless Hamiltonian H ≡ mL 2H / 2 . Here The Schrödinger equation with this Hamiltonian is given by where t = t /(mL 2 ) witht is unscaled time.
In x 1 − x 2 space, the Hamiltonian looks like figure 1. The dotted lines represent delta potentials from trap and interaction between particles. From now on we denote the region where both particles inside the trap as region (1), where one particle in and one particle out of the trap as region (2), and both particles out of trap as region (3).

III. NUMERICAL SIMULATION OF TWO BOSON DECAY
Now, we set up the decay of two interacting identical boson in this Hamiltonian. We choose the initial condition that both particles inside the delta trap as the two boson ground state of V → ∞ case. Specifically, this initial state ψ ini (x 1 , x 2 ) is given by [20] ψ ini (x 1 , with k 1i and k 2i satisfy the equation FIG. 1: Hamiltonian in x 1 − x 2 space. The initial wave function in 0 ≤ x 1 ≤ x 2 ≤ 1 region is obtained from the boson symmetry condition ψ ini (x 1 , x 2 ) = ψ ini (x 2 , x 1 ). In other regions the initial wave function is zero. The normalization constant N ini is chosen to satisfy dx 1 dx 2 |ψ ini | 2 = 1. k 1i and k 2i versus interaction strength g is shown in figure 2. For zero interaction both k 1i and k 2i are same as π, the single particle ground state wavevector. For nonzero g they rapidly deviate form π as g increases, and k 1i approaches to 2π and k 2i approaches to π ( Figure 2) . approaches to 2π and k 2i approaches to π.
To analyze the decay of interacting bosons, we solve the Schrödinger equation directly.
The Schrödinger equation and its formal solution are We use Crank-Nicolson method to solve this equation numerically [21,22].
For the numerical representation of Hamiltonian, we choose the triangular region 0 ≤ , with X max is large enough that in our observing time very little decay products reach near X max . This region is discretized by dividing X max by N x , and all points in the triangular region are arranged in one column vector. The Hamiltonian matrix obtained by discretization of x space and using a finite-difference formula for the second derivatives can be quite large, but it is a sparse matrix as most elements are zero.
For small dt, This method is second order in dt and unitary (i.e. probability is conserved). This is an implicit method, since it contains the inverse operator. The matrix inversion is efficiently implemented by solving linear equations. The time evolution of the wave function is obtained by iterating equation (18).
For the simulations in next sections, the following parameters are used. X max = 45, The convergence of the numerical solutions is checked by changing spatial grid size and time step. We also check numerical simulation with known analytic solutions for special cases g = 0 and g = ∞. To see the effect from the reflection of waves at the boundary the results are examined by changing X max and by putting absorbing potentials near X max . In our parameter regime, those effects are very small and do not change the main results.

IV. RESULTS AND ANALYSIS
For the understanding of the decay of two interacting bosons, a good starting point is the parameter region where we know the exact analytic solutions. In our case, we know exact eigenfunctions of Hamiltonian for two extreme cases, g = 0 and g = ∞. In those cases, the two particle eigenfunctions are obtained by the combination of one particle eigenfunctions, which are known in analytic form. For arbitrary g > 0, the results lie between these two extremes, and the exact analytic forms are not known.

A. Vanishing and infinite interaction limits
When g = 0 there is no interaction between two particles. They act independently, with only a symmetric wavefunction condition. The eigenfunction is written as where the total eigenenergy is E = (k 2 1 + k 2 2 )/2 and φ(k, x) is the one particle eigenfunction with eigen wavevector k. In our model the explicit form of φ is given by where The one particle decay rate can be calculated by outgoing boundary condition, setting the coefficient of outgoing wave c 3 (k) = 0 and solving for k (this is also the pole of scattering matrix). The equation c 3 (k) = 0 has complex solutions, each of them corresponds to different decay modes. We denote the complex solution of c 3 (k) = 0 as k z0 , k z1 , ... with k z0 the lowest decay mode and k z1 next lowest decay mode, etc. For the V = ∞ ground state initial the dominant decay mode is k z0 . Since the decay mode wavefunction is also a complex eigenfunction, its time dependence is given by e −iEt , where E = k 2 z0 /2. The time evolution of one particle probability inside the potential trap P 1in (t) follows the exponential decay where k z0r and k z0i are the real and imaginary parts of k z0 , respectively.
The decay of two interacting bosons in the special cases of g = 0 and g = ∞ is obtained from the single particle decay patterns, respectively.
For g = 0, the two particle wave function is the product of one-particle wave functions, and their decay is just the product of the individual decay. With the condition that the initial wave function was the ground state of V = ∞: If we write the probability of both particle inside the trap as P 1 , probability of one particle in and one out as P 2 and both particles out as P 3 , their dominant time evolutions are Another case we know the exact eigenfunction of Hamiltonian is g = ∞ case. In this case, the two particle eigenfunction is written as and ψ(k 1 , k 2 , x 1 , x 2 ) = ψ(k 1 , k 2 , x 2 , x 1 ) for the x 1 < x 2 region. Like in the case of fermions the probability density is zero along x 1 = x 2 line. The V 0 = ∞ ground state initial condition is given by and ψ ini (x 2 , x 1 ) = ψ ini (x 1 , x 2 ) for x 1 < x 2 . The decay of g = ∞ two bosons involves two decay mode, with lowest wavevector k z0 and next lowest one k z1 . The time evolutions of P 1 , P 2 and P 3 are with and k zjr and k zji are real and imaginary parts of k zj , respectively.

B.
Arbitrary g > 0 case For the general case of 0 < g < ∞ exact analytic eigenfunctions are not known. We use the numerical method of section III to obtain the decay of probabilities. To conduct the simulation, first the initial condition was chosen as the ground state of trap potential strength V = ∞ limit.
Quite interestingly, the numerical results in this section show that a rather simple model can be used to explain interacting boson decay. For the decay of interacting particles, it is expected that the number density of particles shows non-exponential decay. When there are more particles inside the trap it decays faster, and with less particles the decay is slower.
But if we examine the probability P 1 of both particles inside and the probability P 2 of one particle inside and another out separately, they show quite distinctive features.
If we plot the logarithm ln P 1 vs time for various interaction strength g, the graphs show straight lines, meaning the decay is exponential. Furthermore, the decay rate can be obtained by theoretical estimation. Like the decay rate calculation of one-particle case, we can apply the outgoing boundary condition for the wavefunction in region (1). Since the probability of both particles escaping simultaneously is very small due to the repulsive interaction, it is ignored. Then the outgoing boundary condition from region (1) to region (2) can be written as follows.
First, the wavefunction in region (1) ψ (1) , satisfying the Bethe ansatz and the boundary conditions at x 1 = 0 and x 1 = x 2 , can be written as [the form of the coefficients without normalization is given in equations (9) to (12)] and the outgoing eigenfunction in region (2), ψ (2) , can be written as with the boundary condition Conditions (39) and (40) yield four equations with four unknowns B 1 , B 2 , k 1 and k 2 . Solving for k 1 and k 2 we get two equations and two complex wavevectors k 1g and k 2g for their solutions. When we write real and imaginary parts of complex eigenvectors as k 1g = k 1gr + ik 1gi and k 2g = k 2gr + ik 2gi , both of their imaginary parts are negative. Considering that the time evolution of an energy eigenfunction follows e −iEt like the one particle decay mode, it can be expected that exp(−i(k 2 1g + k 2 2g )t/2) dominates in time evolution. When we compare the probability of both particle inside indeed we see that this is what happens. Both are in very good agreements as shown in figure 3. P 1 (t) decays exponentially with the decay rate predicted by outgoing boundary conditions.
(44) Figure 4 shows γ g change for various g. γ g changes a lot for small g, and approaches to the decay rate of g = ∞ slowly. Figure 3 shows the comparison between −γ g t line from theoretical prediction and log P 1 from numerical simulation. They match very well well for all g > 0, thus showing P 1 decays exponentially even with interaction between bosons. Next we consider the time evolution of P 2 , one particle in and one particle out of trap probability. It is more complicated than that of P 1 , since it contains probability inflow from region (1) and outflow into region (3). Like P 1 case, we already know the dominant parts of P 2 (t) for special cases, g = 0 and g = ∞.
For g = 0, the decay of P 2 (t) has the form and for g = ∞, where γ z0 , γ z1 are the lowest and next lowest decay rate of one particle in the potential trap.
For the g = 0 case, both bosons decay from the same mode independently. For the g = ∞ case, two bosons decay from the separate decay modes without interfering since they are almost orthogonal. For general 0 < g < ∞, the time evolution of P 2 (t) will be between (45) and (46) and as g is increased P 2 (t) will change from (45) to (46). We try different models for two regimes where g is not large (weak or moderate repulsion) and where g is very large (strong repulsion), and investigate regions of validity for each model.
For the weak or moderate repulsive interaction, we try a simple model for P 2 decay. If we assume that the probability of both particle escaping simultaneously is very small, which is satisfied when the decay rate is small and interparticle interaction is repulsive, then the inflow from region (1) has very simple form since the dominant part of P 1 satisfies (43) and almost all escaping probability from region (1) goes to region (2). We can write P 2 as where F in is the probability inflow from region (1) to region (2) and F out is the probability outflow from region (2) to region (3). represent P 1 decay rates of g = 0 and g = ∞ cases, respectively. F in is simply γ g e −γg t , which is P 1 (t) outflow from region (1). For the form of outflow F out , we try exponential decay model. In that case, F out is set as −γ 23 P 2 where γ 23 is the decay constant from region (2) to region (3). With this assumption, the solution of (47) has the form P 2 (t) = γ g γ g − γ 23 (e −γ 23 t − e −γgt ).
The decay constant γ 23 is yet undetermined, so (48) becomes one parameter fitting model. This exponential decay model of F out implies that remaining particle in the trap will decay exponentially like one particle decay after other one escapes, with only one decay mode.
Compared with numerical simulation, model (48) shows very good agreements. Furthermore, it shows that even for larger g the fitted parameter γ 23 is very close to the lowest single particle decay rate γ z0 . Figure  for g = 10. k 1i and k 1g are closer to second decay modes, but still P 2 decay to region (3) is dominated by single particle lowest decay rate. Figure 6 shows γ 23 compared to γ z0 and their relative differences for various g. It shows that the relative difference between γ 23 and γ z0 are less than 2% for 0 < g < 17, and the difference increases and approaches to 10% for larger g.
In the strongly repulsive interaction region where the difference between γ 23 and γ z0 increases, the deviation of model 48 from numerical simulation also increases. In this region we try different model which is close to (46). Physical meaning of (46) is that there are two decay modes and two decay modes decay independently without interfering each other. In our case, we have two complex eigenvector k 1g and k 2g from (41) and (42). Assuming that P 2 decays from each complex wavevector and each mode do not interfere each other, we write decay model of P 2 (t) for large g as where This model works better for larger g than (48) as figure 7 shows. The (49) model describes the peak of P 2 (t) well, and discrepancy with the numerical data becomes smaller for larger g.
To see the agreements between numerical simulation and fitting model quantitatively, we consider the absolute mean of relative error η between the probability calculated by numerical simulation P 2num and the probability calculated by model P 2model , defined as with t i s are taken from t = 0.1 to t = 5 by 0.1 intervals. Figure 8 shows plots of η versus interaction strength g plots for two different models.
Model (48) shows good agreements with simulation up to g = 17, with relative error less than 3%. The error increases steadily, reaching more than 5% when g > 60 (seen from figure 2, this is strongly repulsive region). Second model (49) shows large deviation from the numerical simulation for smaller g, but agreements with the simulation becomes better than that of (48) model for g > 67.
Finally, The probability of both particles outside, P 3 , is easily calculated since P 1 + P 2 + P 3 = 1. So total decay mechanism can be described by (43), (48) or (49).     but for larger g the decay rate changes from faster to slower ones.
If we calculate the number density N in (t) inside the potential, it is written as and in our case it simply becomes N in (t) = 2P 1 (t) + P 2 (t). Figure 9 Shows the logarithm of N in versus time. For g = 0, 1 the decay of N in (t) is close to exponential (ln N in close to straight line) but for larger g it is more visible that the decay rate changes from faster to slower ones, as expected.

V. CONCLUSION
We have calculated the decay of two repulsively interacting bosons, initially in the ground state of a potential trap, by numerical simulation. We have found an exponential decay mode for the probability of both bosons inside the trap and have estimated its decay rate theoretically. By applying outgoing boundary condition for the loss of single particle from the trap, we obtain two complex wavevectors corresponding to the two particles inside the trap and corresponding decay rate. The agreement between numerical simulation and theoretical estimation in time evolution of decay probabilities is very good. For describing the probability to have a single particle inside and one particle outside, two simple models were proposed. For small and moderate g, we apply a model in which the remaining particle decays exponentially, whereas for larger g (strongly repulsive) we propose another model in which the modes of each complex wave vector decay separately. Our numerical simulations show very good agreement for weak and moderate interactions with the first model. For stronger interactions, where fermionization effects become relevant, separate exponential decay model becomes necessary and agrees well with simulations. The number density shows that the decay rate changes over time from fast to slower decay for large g. The results show that simple models describe the overall decay mechanism of repulsively interacting bosons well.