2 τ 3 ( GHZ )=1. It also appears a class of entangled three- qubit states represented by | W 〉 state with τ 3 vanished, i.e., κ is a large constant, to ensure R 2 ∼ 0. N$_{c}$ = 8 is the dimension of three-particle basis set. N$_{p}$ is the number of partition in the decomposition. We can easily see that N$_{p}$ determines the number of variables N$_{v}$ = N$_{p}$ +2 N$_{c}$N$_{p}$ . In Eq.13, κ should be large enough to ensure a small residual R 2 to numerically satisfy the constrain given in Eq.11. However, large κ will make the “energy” surface of E very rough, i.e., E ( { p$_{i}$, Ψ$_{i}$ } ) has many local min- ima separated by large “energy” barriers, which causes great difficulties in finding the global minima of three tangle using traditional constrained simulated annealing (CSA) method, because it tends to be trapped at some local minima. Here, we adopt the replica exchange (also known as parallel tempering [18]) Monte Carlo method [12] which simulates M replicas simultaneously each at a different temperature β 0 = 1 /T max < β 1 < ... < β$_{n}$$_{−}$ 2 < β$_{M}$$_{−}$ 1 = 1 /T min covering a range of interest. Each replica runs independently, except that after cer- tain steps, the configurations can be exchanged between neighboring temperatures, according to the Metropolis criterion, The three-tangle of mixed state can be obtained via convex-roof extension[1, 5, 17]. Suppose the density ma- trix of a mixed state ρ can be decomposed into sum of some pure states π$_{i}$ , i.e., where, | Ψ$_{i}$ 〉 is the wave function of a normalized pure state. The three-tangle of the mixed state is defined as the average pure-state concurrence minimized over all possible de- compositions. It is very difficult to develop an universal analytical method to calculate the three-tangle of mixed three- qubit state. So far, there are very limited examples of mixed states whose three-tangle have been obtained analytically[8–11]. Alternatively, one could resort to the numerical methods. To get the three-tangle of mixed states is a constrained minimization problem. The pure state wavefunction can be expanded on the orthonormal three-particle basis, i.e., | Ψ$_{i}$ 〉 = ∑ α c$_{iα}$ | Φ$_{α}$ 〉 . The first constrain is, where Δ H = − ( β$_{i}$ − β$_{i}$$_{−}$ 1 )( E$_{i}$ − E$_{i}$$_{−}$ 1 ), in which E$_{i}$ and E$_{i}$$_{−}$ 1 are the energy of the i -th and i -1-th replica. During the exchange, the detailed balance in canonical ensemble is satisfied. Importantly, the inclusion of high- T config- urations ensures that the lower- T systems can access a broad phase space and avoid becoming trapped at local minima. The replica temperatures are adjust so that the exchange rate between the replicas are all about 20%[19]. while keep the highest temperature β 0 = 1 /T max and lowest temperature β$_{M}$$_{−}$ 1 = 1 /T min fixed. This can be done by adaptively adjusting the temperature of each replica using a recursion method proposed by Berg[20]. We slightly revised the the algorithms to give better per- formance. Suppose, in the n -th iteration, the temper- ature of the i -th replica is β n i and the acceptance rate between β n i and β n i$_{−}$ 1 is a r,n i . In the n +1-th iteration, β n +1 i can be updated as, where, where ρ$_{αβ}$ is the density matrix in the basis of {| Φ$_{α}$ 〉} . Furthermore, the coefficients have to satisfy the normal- ization conditions, The above constrains can be enforced in the simulation by explicitly applying normalization factors at each MC step, whereas the constrain Eq. 11 can be enforced via a penalty function. The “energy function” (target func- tion) reads, c is an empirical parameter that can be adjusted to ac- celerate the convergency. In the present case, we found that by choosing proper c < 1 (in this case c=0.7) the recursion scheme get a suitable β distribution within tens of iterations. We then fixed the replica temperatures for the simulations. Based on convex-roof extension, Lohmayer et al. have provided a complete analysis of mixed three-qubit states composed of a GHZ state and a W state, obtaining where { p$_{i}$, Ψ$_{i}$ } realizing the mixed state as given in Eq.9. R 2 is the residual between the searched density matrix and the target density matrix ρ , defined as | W 〉 = 1 √$_{3}$ ( | 100 〉 + | 010 〉 + | 001 〉 ) (7) ρ = ∑ i p$_{i}$π$_{i}$ , (8) π$_{i}$ = | Ψ$_{i}$ 〉〈 Ψ$_{i}$ | . (9) τ 3 ( ρ ) = min ∑ i p$_{i}$τ 3 ( π$_{i}$ ) . (10) w = { 1 Δ H < 0 e $_{−}$Δ H otherwise , (15) ∑ i p$_{i}$ c$_{i,α}$c ∗ i,β = ρ$_{αβ}$. (11) ∑ α | c$_{iα}$ | 2 = 1 , and ∑ i p$_{i}$ = 1 , p$_{i}$ ≥ 0 . (12) β n +1 i = β n +1 i$_{−}$ 1 + (1 − c + ca n i )( β n i − β n i$_{−}$ $_{1}$) , (16) a n i = a r,n i ( β n M$_{−}$ 1 − β n $_{0}$) ∑ M$_{−}$ 1 i =1 a r,n i ( β n i β n i$_{−}$ $_{1}$) , (17) E ( { p$_{i}$, Ψ$_{i}$ } ) = ∑ i p$_{i}$τ 3 (Ψ$_{i}$) + κR $^{2}$, (13) R 2 = N$_{c}$ ∑ α =1 N$_{c}$ ∑ β =1  $^{N$_{p}$}$ ∑ i =1 p$_{i}$ c$_{i,α}$c ∗ i,β − ρ$_{αβ}$  $^{2}$ . (14)