Quantum entanglement in coupled harmonic oscillator systems: from micro to macro

All of our models, except one, are composed of two groups of identical bosonic oscillators with bilinear interaction, with the state of the total system being a Gaussian pure state. That is, the Hamiltonian is

$$\begin{eqnarray}&&H=\displaystyle \sum _{i=1}^{2}\displaystyle \sum _{{j}_{i}=1}^{{N}_{i}}\left[\displaystyle \frac{{p}_{i,{j}_{i}}^{2}}{2{m}_{i}}+\displaystyle \frac{1}{2}{m}_{i}{\omega }_{i}^{2}{x}_{i,{j}_{i}}^{2}\right]+\alpha \left(\displaystyle \sum _{{j}_{1}=1}^{{N}_{1}}\displaystyle \sum _{{j}_{2}=1}^{{N}_{2}}{x}_{1,{j}_{1}}{x}_{2,{j}_{2}}\right)+\displaystyle \sum _{i=1}^{2}\displaystyle \sum _{\underset{{j}_{i}\gt {l}_{i}}{{j}_{i},{l}_{i}}}^{{N}_{i}}{\beta }_{i}{x}_{i,{j}_{i}}{x}_{i,{l}_{i}},\end{eqnarray} \tag{ 1 }$$

where ${x}_{i,{j}_{i}}$ represent the displacement of the jth particle in the i-th group. ${p}_{i,{j}_{i}}$ be the conjugate momentum of ${x}_{i,{j}_{i}}.\alpha$ and ${\beta }_{i}$ represent the coupling strengths of oscillators between groups 1 and 2 and within group i respectively. The schematic graph of the model is given in figure 1. The difference between these models is the particle number of each group.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Unless mentioned, the initial states considered in this paper are the ground state for the uncoupled oscillators, i.e.

$$\begin{eqnarray}&&\psi ({\boldsymbol{x}},0)={ \mathcal N }\mathrm{exp}\left[-\displaystyle \sum _{i=1}^{{N}_{1}}\displaystyle \frac{{m}_{1}{\omega }_{1}}{2{\rm{\hslash }}}{x}_{1,i}^{2}-\displaystyle \sum _{j=1}^{{N}_{2}}\displaystyle \frac{{m}_{2}{\omega }_{2}}{2{\rm{\hslash }}}{x}_{2,j}^{2}\right],\end{eqnarray} \tag{ 2 }$$

which is Gaussian and isn't the true ground state of this Hamiltonian. There are occasions where we don't choose states of this form, but they're still set to be of the form of generalized Gaussian state. Because the Hamiltonian is quadratic and the interactions are bilinear among micro degrees of freedom, for any given Gaussian-type initial state, the final state of the system at later time can be determined analytically. Some details can be found in appendix A. Knowing the state at any time t, we can discuss the entanglement dynamics of the system.

It should be emphasized that we only discuss bipartite entanglement between two oscillators, and we adopt the unit ${\rm{\hslash }}=1$ when computing.

Here we introduce two terminologies which will be used later:

• Intra-group entanglement is the entanglement between any two particles of the same group, i.e., entanglement 11 or 22, the meaning of which is explained in figure 1.
• When we choose a particle from each group, the entanglement between these two oscillators is the inter-group entanglement, i.e., entanglement 12.

Next, we introduce the separability criterion and the entanglement measure we will use in this paper.

For a state of a continuous system to be separable, Peres-Horodecki-Simon criterion [34–36] is a necessary condition,

$$\begin{eqnarray}&&\tilde{{\rm{V}}}+\displaystyle \frac{i{\rm{\hslash }}}{2}{\rm{\Omega }}\geqslant 0,\end{eqnarray} \tag{ 3 }$$

where $\tilde{{\rm{V}}}={\rm{diag}}(1,1,1,-1)\cdot {\rm{V}}\cdot {\rm{diag}}(1,1,1,-1)$ and ${\rm{V}}$ is the covariance matrix, the computation of which can be simplified with some tricks as shown in appendix B. For a Gaussian state, this is a necessary and sufficient condition [36].

By Sylvester's criterion² , a matrix is positive semi-definite if and only if all of its leading principal minors are non-negative. That is, a matrix is positive semi-definite if and only if the determinants of all of its upper-left corner square submatrices and that of the matrix itself are all non-negative. Because ${\rm{V}}+i{\rm{\hslash }}{\rm{\Omega }}/2$ and $\tilde{{\rm{V}}}+i{\rm{\hslash }}{\rm{\Omega }}/2$ only differ in the fourth row and fourth column, and because ${\rm{V}}+i{\rm{\hslash }}{\rm{\Omega }}/2$ is positive semi-definite³ , (3) is equivalent to

$$\begin{eqnarray}&&{\mathsf{E}}\equiv -\mathrm{det}\left(\tilde{{\rm{V}}}+\displaystyle \frac{i{\rm{\hslash }}}{2}{\rm{\Omega }}\right)\leqslant 0.\end{eqnarray} \tag{ 4 }$$

Later on we will just call this determinant ${\mathsf{E}}.$

Logarithmic negativity [38] is defined as⁴

$$\begin{eqnarray}&&{\mathsf{L}}\equiv \mathrm{ln}| | {\rho }^{{T}_{A}}| | ,\end{eqnarray} \tag{ 5 }$$

where ${\rho }^{{T}_{A}}$ is the partial transpose of ρ and $| | B| | \equiv {\rm{Tr}}\sqrt{{B}^{\dagger }B}.$ It was later proven to be an entanglement monotone, which does not increase under general positive partial transpose preserving (PPT) operations [39]. If the partially transposed state is still a physical state, ${\mathsf{L}}=0$ because the density operator is a semi-positive operator whose trace is 1, which implies that a separable state has ${\mathsf{L}}=0.{\mathsf{L}}$ of a Gaussian state can be computed with the covariance matrix [38], and it has been proven that the larger symplectic value doesn't contribute to the entanglement measure [36, 40].

The Peres-Horodecki-Simon criterion implies that for a Gaussian states, it's separable if and only if the covariance matrix of its partial transpose is physically realizable [36, 41], so is the partial transpose itself. Thus, a Gaussian state is separable if and only if both ${\mathsf{E}}$ and ${\mathsf{L}}$ are equal to zero.

When we're interested in the degree of entanglement, we will evaluate its ${\mathsf{L}}$ because it's a genuine entanglement measure. As to ${\mathsf{E}}$, we keep it as a comparison. Even though being only a separability criterion, not a rigorous entanglement measure, ${\mathsf{E}}$ somehow can measure the entanglement of a state to some extent. Its dynamical behavior is somewhat like ${\mathsf{L}}$ as shown in figures displaying entanglement dynamics.

In this model ${N}_{1}={N}_{2}=1,$ essentially a pair of coupled oscillators, which shows some basic features that can be helpful in understanding how the entanglement dynamics of coupled oscillators works.

Because of the existence of coupling, the entanglement changes with time instead of staying at where it is initially. Thus, even though the state is separable at the beginning, it becomes entangled immediately.

Due to periodicity, they can become separable again after a finite time. Thus, entanglement sudden death occurs in this model. The wave function is of the form

$$\begin{eqnarray}&&\psi ({x}_{1},{x}_{2};t)={ \mathcal N }(t)\mathrm{exp}[-{a}_{1}(t){x}_{1}^{2}-{a}_{2}(t){x}_{2}^{2}+{a}_{12}(t){x}_{1}{x}_{2}].\end{eqnarray} \tag{ 6 }$$

From (6), we can see that what governs whether these two particles are separable is ${a}_{12}(t)$. Applying the propagator (32) over an infinitesimal time dt on a wave function whose ${a}_{12}(t)=0,{a}_{12}(t+{dt})$ can not vanish, so ${a}_{12}(t)=0$ only occurs at discrete values of t. The interaction will make them entangled right after they're separable; in other words, disentanglement can not last for a finite interval in this model.

3.1.1. Entanglement strength

Next we examine how each parameter affects the degree of entanglement, for which an effective way of description is needed. The time evolution of ${\mathsf{L}}$ itself is too complicated to be analyzed directly, so instead we take a numerical approach by finding the maximum value of ${\mathsf{L}}$ in a given time interval. Plotting the maximum value of ${\mathsf{L}}$ against a given parameter, this kind of figures can, to some extent, represents how a parameter affects the strength of entanglement. Later on whenever we mention entanglement strength, we are referring to the maximum value of ${\mathsf{L}}$ for a given system that can be found during its evolution, denoted ${{\mathsf{L}}}_{\mathrm{max}}$. Please note that entanglement strength as defined here, is by no means a rigorous measure, but of the many aspects of entanglement dynamics it is the one that we focus on in this paper.

From figure 3(a), we find that the larger the coupling strength is, the stronger the maximum entanglement measure is, manifesting the fact that the entanglement is activated by the interaction.

Then, let us see how frequency and mass affect entanglement strength. Figure 3 seems to imply they lower the entanglement strength, but it' is not the whole story. The model has an initial state as (2), which is influenced by frequency and mass as well. It is their roles both in the initial state and the Hamiltonian that causes them to decrease entanglement strength. How they actually influence entanglement strength is superfluous when the degrees of freedom becomes larger and larger as we will see in the following sections.

Zoom In Zoom Out Reset image size

Here ${N}_{1}=1$ and ${N}_{2}=2$. It has just one extra oscillator compared to the one-to-one model, yet it makes a huge difference. We neglect the subscript 2 of ${\beta }_{2}$, denoting β.

3.2.1. Environment and entanglement

3.2.2. Entanglement strength

The first two diagrams of figure 4 show entanglement strength versus interactions. Here we can see drops and rises, unlike in the one-to-one model. The increase of α increases the entanglement strength of ${{\mathsf{L}}}_{12}$. The increase of β increases the entanglement strength of ${{\mathsf{L}}}_{22}$ but decreases ${{\mathsf{L}}}_{12}$. This can be understood as follows:

Zoom In Zoom Out Reset image size

Consider a bipartite system with two parties 1 and 2. Since two separable parties can not be entangled without interaction, unsurprisingly a stronger interaction results in stronger entanglement. There are some parameters that can enhance the entanglement strength, one of which is interaction, while there are some that do the reverse. F us call them entanglement enhancers and reducers.

Now, let us add an environment and treat the whole environment as a party E. Entanglement can be established directly by the interaction between 1 and 2, or indirectly via the environment. When the indirect interaction was turned off, the entanglement strength that would be achievable by the direct interaction alone will be called direct linkage; while if the direct interaction was switched off and the entanglement between 1 and 2 thus would be created via the environment indirectly, the entanglement strength then is called indirect linkage. Entanglement enhancers for 1 and 2 can strengthen the direct linkage, while enhancers for 1 and E, 2 and E, and enhancers for the entanglement within E increase the indirect linkages, similarly for reducers.

Both types of linkages contribute to the entanglement strength between 1 and 2, but their effects don't add up: when one is dominant in a given context, strengthening the other will often decrease the entanglement strength. The graph of ${{\mathsf{L}}}_{\mathrm{max}}$ versus a parameter usually can tell us which one holds the rein. For example, if the indirect linkage dominates at a given set of parameters, then the entanglement strength will drop if we enhance the direct linkage a bit. But since sometimes the trend of change is unclear, we can not always be sure of which is dominant. Note that they play unequal roles in determining entanglement strength. Generally, the direct linkage has more weight in this aspect, and it often results in more explicit change in entanglement strength than the indirect one.

3.2.3. Disentanglement

Figure 5 clearly shows that under certain circumstances, entanglement 22 may disappear for a duration, unlike the one-to-one model, where the disentanglement only appears instantly. Due to the periodicity of the system, this disentanglement can not last forever. We will call it finite-life disentanglement from now on, in contrast to finite-time disentanglement [5], or entanglement sudden death [7], which is defined by entanglement disappearing in a finite time rather than asymptotically. Finite life embodies the facts that the disentanglement doesn't last forever and that this phenomenon happens in a finite time. Thus, it is just a special kind of entanglement sudden death. Note that this kind of phenomenon was already discovered [9–11], but to describe it specifically, introducing a simple term seems better than a lengthy description. On the other hand, finite-life disentanglement doesn't happen between 12, on which we provided an argument in [42].

Zoom In Zoom Out Reset image size

It is interesting to know when finite-life disentanglement can occur between 22. If finite-life disentanglement does happen, ${\mathsf{E}}\leqslant 0$. However, being differentiable to any order, ${\mathsf{E}}$ should be less than zero during this time interval. Thus, we can find the occurrence of finite-life disentanglement by finding if ${\mathsf{E}}$ can be less than zero at some time. Here we draw some graphs, as in figure 6, which shows that if the direct linkage becomes weaker or if the indirect linkage becomes stronger, then finite-life disentanglement will be more likely to occur.

Zoom In Zoom Out Reset image size

The Hamiltonian for this model is

$$\begin{eqnarray}&&H=\displaystyle \sum _{i=1}^{4}\left(\displaystyle \frac{{p}_{i}^{2}}{2m}+\displaystyle \frac{1}{2}m{\omega }^{2}{x}_{i}^{2}\right)+\alpha ({x}_{1}{x}_{2}+{x}_{3}{x}_{4})+\beta ({x}_{1}{x}_{3}+{x}_{2}{x}_{4})+\gamma ({x}_{1}{x}_{4}+{x}_{2}{x}_{3}).\end{eqnarray} \tag{ 7 }$$

figure 7 gives a schematic graph of the boxedcross model. Its interaction pattern, governed by three parameters, makes this model very interesting. Note that in this model a label is given to each oscillator since they can not be separated into two groups now. Again, we set their states as the initial states of decoupled oscillators, and because they have identical mass and frequency, their initial states are the same.

Zoom In Zoom Out Reset image size

Figure 8 shows finite-life disentanglement appearing in this model. figure 9 then exhibits how each interaction parameter changes the entanglement strength, where competition between direct and indirect linkages can be observed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Because the model is more symmetric than the former one, an empirical relation following which finite-life disentanglement happens can be observed, as shown in figure 10. This figure demonstrates that finite-life disentanglement only happens when the summation of the indirect interactions is stronger than the direct interaction. What we mean by direct and indirect interactions, is that, for example, for the pair 1 and 2, its direct interaction is α, by the Hamiltonian term $\alpha {x}_{1}{x}_{2}$, and its indirect ones are β and γ, because 1 and 2 can interact indirectly by the coupling β and γ, as mentioned in section 3.2.1. Thus, for the pair 1 and 2 to have finite-life disentanglement, $\beta +\gamma$ has to be larger than α. Once again, it demonstrates a stronger indirect linkage and a weaker direct linkage make finite-life more likely to happen.

Zoom In Zoom Out Reset image size

The condition that the summation of indirect interactions should be larger than the direct one is only a necessary condition, not a sufficient one. Figure 8 is a good example for illustration. Its three coupling parameters satisfy the triangle inequality, but we only find finite-life disentanglement between 12.

Here we consider the case where there is only one group of identical particles, due to which we drop all the subscripts that specify the groups. We found that, if the particle number N is large enough, it sets an approximate upper bound, $1/N$, for entanglement.

The state of our system is always of this form

$$\begin{eqnarray}&&\psi ={ \mathcal N }\mathrm{exp}\left[-a\displaystyle \sum _{i=1}^{N}{x}_{i}^{2}+2b\displaystyle \sum _{i,j\gt i}^{N}{x}_{i}{x}_{j}\right],\end{eqnarray} \tag{ 8 }$$

where a and b are complex. By assuming a and b are real, we can derive⁶

$$\begin{eqnarray}&&\mathrm{sup}{\mathsf{L}}=\displaystyle \frac{1}{2}\mathrm{ln}\left(1+\displaystyle \frac{2}{N-2}\right)\simeq \displaystyle \frac{1}{N-2}\simeq \displaystyle \frac{1}{N}.\end{eqnarray} \tag{ 9 }$$

Note that $1/(N-2)$ is also an upper bound as shown in appendix C.3. When $b\lt 0$,

$$\begin{eqnarray}&&{\mathsf{L}}\simeq \displaystyle \frac{1}{N-r}{\rm{}}\;{\rm{at}}\;{\rm{large}}\;{\rm{}}N,\end{eqnarray} \tag{ 10 }$$

where $r\equiv a/b$.

What does the above result imply? First, the maximum entanglement achievable by this system is bounded (for $N\ne 2$), which decreases as N increases; second, at larger N, by (10) in a wider range of r will ${\mathsf{L}}$ be near $1/N$, implying in the whole state vector space the subspace where ${\mathsf{L}}\sim 1/N$ becomes larger. Even though a and b are in general complex, the existence of a least upper bound is consistent with the numerical data to be shown below.

In figure 11, the coupling is set extremely tiny. When the initial state is very close to the ground state of decoupled oscillators, i.e. ${a}_{0}=m\omega /(2{\rm{\hslash }})$ and ${b}_{0}=0$, its maximum entanglement measure is constrained by the weak interaction. There is a region of a and b for the initial state where the maximum entanglement stays at the order of the coupling. This means that, by (10), the ratio $r=a/b$ stays large (assuming they're real) if the coupling is weak. As the initial state moves away from this region, its maximum entanglement starts to climb, up to about $1/N$. Interestingly, even when the ratio ${r}_{0}={a}_{0}/{b}_{0}\gg N$, the entanglement strength still grows to $1/N$ if the a₀ is far enough from $m\omega /(2{\rm{\hslash }})$ or b₀ far enough from 0, namely, $r=a/b$ can become much smaller than N as the state evolves.

Zoom In Zoom Out Reset image size

Now let us look at figure 12, where we show how the entanglement strength changes as the interaction β changes. Intuitively, when the interaction becomes larger, we expect it to enhance the entanglement strength. In terms of real a and b, this means that a larger β can yield a smaller $r=a/b.$ Interestingly, from this figure we can observe that the entanglement strength stops growing after it reaches $1/N$, which is exactly the upper bound we found by assuming a and b are real. Note that during the state evolution a and b are rarely real. Therefore, even though a and b are in general complex, the exact upper bound should be very close to or same as the upper bound we proposed, as shown by figure 12(b).

Zoom In Zoom Out Reset image size

The state of this system is

$$\begin{eqnarray}&&\psi ({{\boldsymbol{x}}}_{1},{{\boldsymbol{x}}}_{2})={ \mathcal N }\mathrm{exp}\left[-\displaystyle \sum _{i=1}^{{N}_{1}}{b}_{11}{x}_{1,i}^{2}+2{b}_{1}\displaystyle \sum _{i\gt j}^{{N}_{1}}{x}_{1,i}{x}_{1,j}-\displaystyle \sum _{i=1}^{{N}_{2}}{b}_{22}{x}_{2,i}^{2}+2{b}_{2}\displaystyle \sum _{i\gt j}^{{N}_{2}}{x}_{2,i}{x}_{2,j}+{c}_{12}\displaystyle \sum _{i=1}^{{N}_{1}}\displaystyle \sum _{j=2}^{{N}_{2}}{x}_{1,i}{x}_{2,j}\right],\end{eqnarray} \tag{ 11 }$$

for two groups of bosons. In this part we deal only with intra-group entanglement. In previous section we have shown that for the one-group model $\mathrm{sup}{\mathsf{L}}\simeq 1/N$. Here we find similar behavior, namely, for the intra-group entanglement ${{\mathsf{L}}}_{{ii}}$:

$$\begin{eqnarray}&&\mathrm{sup}{{\mathsf{L}}}_{{ii}}=\displaystyle \frac{1}{2}\mathrm{ln}\left(1+\displaystyle \frac{2}{{N}_{i}-2}\right)\simeq \displaystyle \frac{1}{{N}_{i}-2}\simeq \displaystyle \frac{1}{{N}_{i}},{\rm{}}\;{\rm{for}}\;{\rm{large}}\;{\rm{enough}}\;{\rm{}}{N}_{i},\end{eqnarray} \tag{ 12 }$$

as discussed in appendix C.2 and C.3.

Figure 13(a) shows that ${({{\mathsf{L}}}_{{ii}})}_{\mathrm{max}}$ stops growing after it reaches $1/{N}_{i}$. Figure 13(b) then shows ${{\mathsf{L}}}_{{ii}}$ can not surpass $\mathrm{sup}{{\mathsf{L}}}_{{ii}}$. We can see the difference is very tiny, but ${{\mathsf{L}}}_{{ii}}$ just can not exceed the supremum. Finally, figure 13(c), like before, demonstrates $1/{N}_{i}$ is just an approximation of $\mathrm{sup}{{\mathsf{L}}}_{{ii}},$ so ${{\mathsf{L}}}_{{ii}}$ can be larger than $1/{N}_{i}.$ It is still a good approximation for large N_i nonetheless.

Zoom In Zoom Out Reset image size

In this subsection we will discuss a special case of the two-group model, where group 1, as the system, consists of two oscillators and group 2, as the environment, has a large number of oscillators. We're interested in the entanglement within the system, and that between the system and the environment, i.e. entanglement 11 and 12. The entanglement dynamics of a system consists of two or three oscillators under the influence of a common environment has been studied in work like [43–50]. The phase diagram for general environment, the nonresonant oscillator under weak and strong dissipation, the high-temperature entanglement and the condition for non-equilibrium steady states have been investigated to some extent. Here we would like to add one new perspective — how does the size of the environment affect the entanglement among the system. Unlike previous works, the degrees of freedom of our environment is discrete and finite, so our system never becomes steady, and thus instead of just focusing on the long-time behavior we pick a large number of random time points⁷ to understand its entanglement property. Calculating how much percentage of the time points when the oscillators in the system are disentangled would tell us how frequent the state of the system becomes separable. We limit the discussion on this model to how α and N₂ affect the proportion of time during which the parties under consideration, either 11 or 12, is separable (abbreviated as 'proportion' afterward), shown in figures 14 and 15. Below we will list the trends we observed in the parameter regime we scanned. For entanglement 11

• 1.  
  Increasing both N₂ and α makes the system more prone to become separable.
• 2.  
  No matter how large N₂ is, α has to be large enough to make the system disentangled.
• 3.  
  If N₂ is large enough, at around ${10}^{2}\sim {10}^{3}$ in this case, further increasing the size of the environment has little influence on the proportion. Likewise, at large N₂ and with a large enough α, further enhancing α doesn't make a noticeable difference.
• 4.  
  Because increasing either α or N₂ matters little if each is large enough, we conjecture that at given $({\beta }_{1},{\beta }_{2},{\omega }_{1},{\omega }_{2},{m}_{1},{m}_{2})$ and with a specific initial state, there is a maximum for the proportion in this parameter regime. For the system shown in figure 14 it is about $34 \% .$

Zoom In Zoom Out Reset image size

As for entanglement 12

• 1.  
  At small N₂, the proportion does not have clear tendency. If N₂ is large enough, the proportion increases as N₂ gets larger, and it is almost always separable at very large N₂.
• 2.  
  At small N₂, when changing α we can not really observe an obvious trend. However, at larger N₂, we can clearly see a drop of proportion around some value of α. Then the proportion climbs up again as α becomes even larger. With very large N₂, the proportion is very high regardless of $\alpha .$
• 3.  
  Intuitively, stronger α might have made 12 harder to be separable, which turns out not to be true. It may be because α also makes the indirect linkage stronger. There's a range of value of α that makes a fine balance between direct and indirect linkages. Beyond that value, the indirect linkage dominates.
• 4.  
  As pointed out in section 4, when ${N}_{2}=1$ finite-life disentanglement doesn't happen between 12. However, when N₂ is large, it's hard to find when 12 is entangled, which clearly shows the effect of ${N}_{2}.$

We should emphasize that due to the nature of randomness, at different runs the results are likely to be different, but not significantly. Besides, there's also no strong reason that a larger N₂ or α 'must always' makes 11 or 12 separable more often or not, but the overall tendency can still be found. Please note that in order to keep the eigenfrequencies real and positive, α can not exceed a certain value, which is specific for each system. Therefore it is impossible to take α to infinity without altering other parameters.

We estimate the entanglement strength by taking the maximum value of the entanglement measure, ${{\mathsf{L}}}_{\mathrm{max}}$, in a finite time. As long as a parameter can change the course of how a state evolves, it can change the entanglement dynamics of a system, but they don't have the same capability in changing the entanglement strength.

To explain how the environment, or the third party, influences entanglement strength, we introduce the idea of direct and indirect linkages in section 3.2.2. In general, a parameter that affects the direct linkage, has more weight in changing entanglement strength than that affecting indirect linkage. Among those influencing the indirect linkage, factors that totally pertaining to the environment are even weaker.

We mentioned that when one linkage is dominant, rising the other often decreases the entanglement strength. We should emphasize that sometimes we can see more than one drops and rises, or no obvious drop on the ${{\mathsf{L}}}_{\mathrm{max}}$ versus some parameter graphs. Thus, the direct and indirect linkages work together in a very complicated manner which we don't fully comprehend, or more likely, interpreting it solely by direct and indirect linkages isn't enough.

As discussed in sections 4.1 and 4.2, ${\mathsf{L}}$ may have an inherent upper bound that is determined by the state itself. Hence, no matter how we change the parameters of this system, the entanglement can not go beyond this bound.

In our models, we find that the systems can disentangle for a duration, and turn entangled again, which we call finite-life disentanglement. This is interesting for two reasons: (1) this revival of entanglement doesn't come immediately after the disentanglement. After all, the initial state itself is often set to be separable. It's not surprising at all if it merely becomes separable again just in instants; (2) our system is constantly evolving and doesn't stabilize to a certain state, which would make the disentanglement last forever if this particular state is separable. That finite-life disentanglement can even happen in models as simple as the one-to-two model, section 3.2, suggests that such a phenomenon may be more common than we've thought.

If the direct linkage is weaker than the indirect linkage, then finite-life disentanglement becomes more likely to occur. Note that unlike for entanglement strength, for finite-life disentanglement direct and indirect linkages play on a quite even ground, as shown in section 3.3, where we found that the condition for finite-life disentanglement to occur is simply that the sum of indirect interactions should be larger than the direct interaction; on the contrary, even when this condition is satisfied, increasing the indirect interactions and thus the indirect linkage, the entanglement strength often climbs very slowly, and sometimes even drops.

We should emphasize this is an oversimplified interpretation. For instance, it can not explain why finite life disentanglement doesn't happen between 12 for the one-to-two model, section 3.2. The complete mechanism should be more complicated.

One more interesting thing to note: when the indirect linkage is much stronger than the direct one, on one hand the entanglement strength is enhanced, but on the other hand finite-life disentanglement also can become more and more common. Hence raising entanglement strength between two parties may be accompanied with more frequent finite-life disentanglement; they're not mutually exclusive.

5.3.1. Supremum

We found $\mathrm{sup}{\mathsf{L}}=1/(N-2)$ in section 4.1, and for the two-group model in section 4.2 $\mathrm{sup}{{\mathsf{L}}}_{{ii}}=1/({N}_{i}-2)$. Is this just a coincidence or is there any deeper reason behind? This could be partially answered if we start by observing a spacial class of Gaussian pure states which can be separated into $| \psi ^{\prime} \rangle =| {\psi }_{1}\rangle | {\psi }_{2}\rangle$, where $| {\psi }_{i}\rangle$ is the state for group i and is also Gaussian. For these particular states,

$$\begin{eqnarray}&&\mathrm{sup}{\mathsf{L}}{(| \psi ^{\prime} \rangle )}_{{ii}}=\mathrm{sup}{\mathsf{L}}(| {\psi }_{i}\rangle ).\end{eqnarray} \tag{ 13 }$$

That is, for these states, the supremum of pairwise entanglement of group i, is equal to the supremum pairwise entanglement of a group of bosons whose state is pure and Gaussian. Hence, for a general two-group Gaussian pure state $| \psi \rangle$

$$\begin{eqnarray}&&\mathrm{sup}{\mathsf{L}}{(| \psi \rangle )}_{{ii}}\geqslant \mathrm{sup}{\mathsf{L}}{(| \psi ^{\prime} \rangle )}_{{ii}}=\mathrm{sup}{\mathsf{L}}(| {\psi }_{i}\rangle ),\end{eqnarray} \tag{ 14 }$$

that is, because we've already known the supremum for a special class of states, the supremum of a general state should be equal or larger than it. As it turns out, they are the same, so the possible extra correlation to another group doesn't change the supremum. Does this hold true for an N-group identical oscillator model? This would be a hard problem, but it may be worth further examinations.

These results were derived under the assumption that the state parameters of the wave function are real. However, accompanied with the numerical data, we think even in the more general scenario where those parameters are real, the upper bounds for both cases should be very close to each other. Our results suggests that the entanglement between two parties in a large system will be reduced according to the system size. However, we set a strong constraint on our systems: identicalness. Because the entanglement for a pure bipartite state is unbounded, given a state of this form $\rho (1,2)\otimes \rho (3,4,\ldots )$, the entanglement between 1 and 2 is unbounded. Thus, that the entanglement is bounded is caused by the identicalness. Besides, the total state is a pure Gaussian state. If it's a mixed Gaussian state, or a more general state, the final results will be possibly different. Still, we think our results are nontrivial, and that similar results, i.e. entanglement scaling with the particle number, possibly in a different way, may be found for other models as well. This can be a route to understanding how entanglement behaves in large systems.

The size of our system is totally under our control. Realistically, the size of a macroscopic quantum system requires some methods to measure and determine [22, 24–26]. Entanglement scaling with the size for other kinds of systems has also been researched, such as area laws [51–53].

It's also worth mentioning an entanglement measure, negativity [38], defined as ${\mathsf{N}}=(| | {\rho }^{{T}_{A}}| | -1)/2$, has a similar relation as ${\mathsf{L}}$

$$\begin{eqnarray}&&\mathrm{sup}{\mathsf{N}}=\mathrm{sup}\displaystyle \frac{\mathrm{exp}({\mathsf{L}})-1}{2}\simeq \displaystyle \frac{\mathrm{exp}(1/N)-1}{2}\simeq \displaystyle \frac{1+\tfrac{1}{N}-1}{2}=\displaystyle \frac{1}{2N}{\rm{}}\;{\rm{at}}\;{\rm{large}}\;{\rm{}}N.\end{eqnarray} \tag{ 15 }$$

5.3.2. Unboundedness when ${N}_{i}=2$

5.3.3. Two to N₂

In section 4.3 we discussed how an environment with a scalable size influences the entanglement within the system and that between the system and environment. Please note the observations are only applicable to the parameter subspace we've investigated.

We found that a larger N₂, thus a larger environment, as well as a larger inter-group interaction α makes 11 disentangled more often. This is because when α and N₂ enhance the indirect linkage, the disentanglement becomes more frequent. We also learned that if the environment is big enough, making it bigger doesn't make much difference to how often the system disentangle; similarly, with a large environment if α is strong enough, making it even stronger has little effect on the frequency of disentanglement.

For entanglement 12, the most interesting finding is at large N₂ it becomes disentangled very often. If N₂ is very large, then it's very hard to find when 12 is entangled. At large N₂, we can also see the competition between direct and indirect linkages. Both direct and indirect linkages grow as α is strengthened, but direct linkages loses when α is above a certain strength.

It would be interesting to compare our model with the quantum Brownian motion (QBM) model for two oscillators in the same environment [43, 54]. In QBM models the environment is a collection of independent oscillators with prescribed oscillation frequencies and spectral density. The dynamics of the reduced density matrix is governed by master equations or Langevin equations, where the environment has already been integrated out.

In our model, there is permutation symmetry among each group and the interaction pattern is different from those of QBM. Due to this symmetry, after diagonalizing group 2 there are ${N}_{2}-1$ degenerate degrees of freedom that are irrelevant to the dynamics of and the entanglement within the system (group 1), as one can see in appendix A. From this perspective enlarging N₂ doesn't actually increase the effective degrees of freedom coupled to the system, but merely alters some parameters. However, in our model we can study not only the entanglement between the system degrees of freedom but also the entanglement between the system and environment as shown in figure 15. Our model may capture some features of entanglement within a system that is coupled to a large but finite environment that possesses distinct degrees of freedom. In addition, entanglements 12 and 22 depend on those degenerate degrees of freedom of the environment, so they are nontrivial when discussing those entanglements.

Zoom In Zoom Out Reset image size

The Hamiltonian is:

$$\begin{eqnarray}&&H=\displaystyle \sum _{i}^{N}\displaystyle \sum _{{j}_{i}}^{{N}_{i}}\left[\displaystyle \frac{{p}_{i,{j}_{i}}^{2}}{2{m}_{i}}+\displaystyle \frac{1}{2}{m}_{i}{\omega }_{i}^{2}{x}_{i,{j}_{i}}^{2}\right]+\displaystyle \sum _{i,k\gt i}^{N}{\alpha }_{i,k}\left(\displaystyle \sum _{{j}_{i}}^{{N}_{i}}\displaystyle \sum _{{j}_{k}}^{{N}_{k}}{x}_{i,{j}_{i}}{x}_{k,{j}_{k}}\right)+\displaystyle \sum _{i}^{N}\displaystyle \sum _{{j}_{i},{l}_{i}\gt {j}_{i}}^{{N}_{i}}{\beta }_{i}{x}_{i,{j}_{i}}{x}_{i,{l}_{i}},\end{eqnarray} \tag{ 16 }$$

which can be written more elegantly as

$$\begin{eqnarray}&&H=\displaystyle \sum _{i}^{N}\left(\displaystyle \frac{1}{2{m}_{i}}{{\boldsymbol{p}}}_{i}^{{\rm{T}}}{{\boldsymbol{p}}}_{i}+\displaystyle \frac{1}{2}{{\boldsymbol{x}}}_{i}^{{\rm{T}}}{{\rm{U}}}_{i}{{\boldsymbol{x}}}_{i}\right)+\displaystyle \frac{1}{2}{{\boldsymbol{X}}}^{{\rm{T}}}{\rm{A}}{\boldsymbol{X}},\end{eqnarray} \tag{ 17 }$$

where ${{\boldsymbol{x}}}_{i}$ comprise coordinates of the i-th group, similarly for ${\boldsymbol{p}}$, and ${\boldsymbol{X}}$ has its elements⁸ ${X}_{i}={({N}_{i})}^{-1/2}{\sum }_{{j}_{i}=1}^{{N}_{i}}{x}_{i,{j}_{i}}$. Also,

$$\begin{eqnarray}{{\rm{A}}}_{i,j}=\left\{\begin{array}{ll}0 & {\rm{if}}\;{\rm{}}i=j\\ \sqrt{{N}_{i}{N}_{j}}{\alpha }_{i,j} & {\rm{otherwise}}\end{array}\right.\end{eqnarray} \tag{ 18 }$$

and

$$\begin{eqnarray}&&{({{\rm{U}}}_{i})}_{{jk}}={\beta }_{i}+{\delta }_{{jk}}({m}_{i}{\omega }_{i}^{2}-{\beta }_{i}).\end{eqnarray} \tag{ 19 }$$

Because the matrices ${{\rm{U}}}_{i}$ are symmetric (hermitian), they certainly can be diagonalized. By the fact that a rank-N matrix ${\rm{M}}$

$$\begin{eqnarray}&&{{\rm{M}}}_{{ij}}=b+{\delta }_{{ij}}(a-b)\end{eqnarray} \tag{ 20 }$$

has eigenvalues a − b and $a+(N-1)b$, the first of them having $(N-1)$-fold degeneracy. The eigenvector corresponds to the eigenvalue $a+(N-1)b$ is of the form $(1,1,\ldots ,1),$ which is exactly the center of mass coordinate multiplied by N.

Due to the degeneracy, we can choose any vectors that are mutually orthogonal and orthogonal to $(1,1,\ldots ,1).$ An easy way to finding them all is using the Gram-Schmidt process. For example, we can start from $\tfrac{1}{\sqrt{{N}_{i}}}(1,1,\ldots ,1),(1,0,0,\ldots ),(0,1,0,0,\ldots )$ , ..., $(0,\ldots ,0,1,0).$ Having all vectors normalized will make our work afterward a bit neater. After all this, we have the transformation matrix:

$$\begin{eqnarray}&&{{\rm{O}}}_{i}=\left(\begin{array}{c}{\hat{v}}_{1}^{{\rm{T}}}\\ {\hat{v}}_{2}^{{\rm{T}}}\\ \vdots \\ {\hat{v}}_{{N}_{i}}^{{\rm{T}}}\end{array}\right),\end{eqnarray} \tag{ 21 }$$

where ${\hat{v}}_{1}^{{\rm{T}}}=\tfrac{1}{\sqrt{{N}_{i}}}(1,1,\ldots ,1)$. This is an orthogonal matrix, and its inverse matrix thus is its own transpose.

Now let us throw all stuff into the pot. From (17), making the transformation ${{\boldsymbol{u}}}_{i}={{\rm{O}}}_{i}{{\boldsymbol{x}}}_{i}$ where ${\rm{O}}$ is defined as in (21) along with the corresponding momenta ${{\boldsymbol{P}}}_{i}={{\rm{O}}}_{i}{{\boldsymbol{p}}}_{i}$:

$$\begin{eqnarray}&&H=\displaystyle \sum _{i}^{N}\left(\displaystyle \frac{1}{2{m}_{i}}{{\boldsymbol{P}}}_{i}^{{\rm{T}}}{{\rm{O}}}_{i}{{\rm{O}}}_{i}^{-1}{{\boldsymbol{P}}}_{i}+\displaystyle \frac{1}{2}{{\boldsymbol{u}}}_{i}^{{\rm{T}}}{{\rm{O}}}_{i}{{\rm{U}}}_{i}{{\rm{O}}}_{i}^{-1}{{\boldsymbol{u}}}_{i}\right)+\displaystyle \frac{1}{2}{{\boldsymbol{X}}}^{{\rm{T}}}{\rm{A}}{\boldsymbol{X}}\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&=\displaystyle \sum _{i}^{N}\left(\displaystyle \frac{1}{2{m}_{i}}{{\boldsymbol{P}}}_{i}^{{\rm{T}}}{{\boldsymbol{P}}}_{i}+\displaystyle \frac{1}{2}{{\boldsymbol{u}}}_{i}^{{\rm{T}}}{{\rm{D}}}_{i}{{\boldsymbol{u}}}_{i}\right)+\displaystyle \frac{1}{2}{{\boldsymbol{X}}}^{{\rm{T}}}{\rm{A}}{\boldsymbol{X}},\end{eqnarray} \tag{ 23 }$$

where ${{\rm{D}}}_{i}$ is

$$\begin{eqnarray}&&{{\rm{D}}}_{i}={\rm{diag}}({m}_{i}{\omega }_{i}^{2}+({N}_{i}-1){\beta }_{i},{m}_{i}{\omega }_{i}^{2}-{\beta }_{i},{m}_{i}{\omega }_{i}^{2}-{\beta }_{i},\ldots ).\end{eqnarray} \tag{ 24 }$$

We can alternatively write the Hamiltonian as

$$\begin{eqnarray}H & = & \displaystyle \frac{1}{2}\displaystyle \sum _{i}^{N}\left[\displaystyle \sum _{j=1}^{{N}_{i}}\displaystyle \frac{{P}_{i,j}^{2}}{{m}_{i}}+\displaystyle \sum _{j=2}^{{N}_{i}}({m}_{i}{\omega }_{i}^{2}-{\beta }_{i}){u}_{i,j}^{2}+({m}_{i}{\omega }_{i}^{2}+({N}_{i}-1){\beta }_{i}){X}_{i}^{2}\right]\\ & & +\displaystyle \sum _{i=1,j\gt i}^{N}\sqrt{{N}_{i}{N}_{j}}{\alpha }_{i,j}{X}_{i}{X}_{j}.\end{eqnarray} \tag{ 25 }$$

Note that ${({u}_{i})}_{1}\equiv {u}_{i,1}={X}_{i}$. In this Hamiltonian, only virtual oscillators X_i have interactions among themselves, while all other oscillators are isolated. Also, all ${u}_{i,j}$ of the same i with $j\gt 2$ undergo the same dynamic. We've successfully reduced this Hamiltonian into one with a much simpler interaction pattern.

Under the same Hamiltonian (17) or (25), we want to solve the wave function at time t, whose initial state is⁹

$$\begin{eqnarray}&&\psi ({{\boldsymbol{x}}}_{1},{{\boldsymbol{x}}}_{2},\ldots ,{{\boldsymbol{x}}}_{N};0)={ \mathcal N }(0)\mathrm{exp}\left[\displaystyle \sum _{i=1}^{N}{{\boldsymbol{x}}}_{i}^{{\rm{T}}}{{\rm{B}}}_{i}{{\boldsymbol{x}}}_{i}+{{\boldsymbol{X}}}^{{\rm{T}}}{\rm{C}}{\boldsymbol{X}}\right],\end{eqnarray} \tag{ 26 }$$

where

$$\begin{eqnarray}&&{({{\rm{B}}}_{i})}_{{jk}}={b}_{i}-{\delta }_{{jk}}({b}_{{ii}}+{b}_{i})\end{eqnarray} \tag{ 27 }$$

and

$$\begin{eqnarray}{{\rm{C}}}_{{ij}}={{\rm{C}}}_{{ji}}=\left\{\begin{array}{ll}0 & {\rm{if}}\;{\rm{}}i=j\\ \displaystyle \frac{1}{2}{c}_{{ij}} & {\rm{otherwise}}\end{array}\right..\end{eqnarray} \tag{ 28 }$$

Making the same transformation as before, we have¹⁰ :

$$\begin{eqnarray}&&{{\boldsymbol{x}}}_{i}\longrightarrow {{\boldsymbol{u}}}_{i}\;:\psi ({{\boldsymbol{u}}}_{1},{{\boldsymbol{u}}}_{2},\ldots ,{{\boldsymbol{u}}}_{N};0)={ \mathcal N }(0)\mathrm{exp}\left[\displaystyle \sum _{i=1}^{N}{{\boldsymbol{u}}}_{i}^{{\rm{T}}}{{\rm{F}}}_{i}{\boldsymbol{u}}+{{\boldsymbol{X}}}^{{\rm{T}}}{\rm{C}}{\boldsymbol{X}}\right],\end{eqnarray} \tag{ 29 }$$

where

$$\begin{eqnarray}&&{{\rm{F}}}_{i}={\rm{diag}}(-{b}_{{ii}}+({N}_{i}-1){b}_{i},-{b}_{{ii}}-{b}_{i},-{b}_{{ii}}-{b}_{i},\ldots ).\end{eqnarray} \tag{ 30 }$$

We can write it in a more lengthy but enlightening way:

$$\begin{eqnarray}\psi (0) & = & { \mathcal N }(0)\mathrm{exp}\left\{-\displaystyle \sum _{i=1}^{N}[{b}_{{ii}}-({N}_{i}-1){b}_{i}]{X}_{i}^{2}\right.\\ & & \left.+\displaystyle \sum _{i=1,j\gt i}^{N}{c}_{{ij}}{X}_{i}{X}_{j}\right\}\mathrm{exp}\left[-\displaystyle \sum _{i=1}^{N}\displaystyle \sum _{j=2}^{{N}_{i}}({b}_{{ii}}+{b}_{i}){u}_{i,j}^{2}\right].\end{eqnarray} \tag{ 31 }$$

The first and second part evolve on their own, according to the Hamiltonian (25). Therefore, this problem practically becomes solving the wave function of N coupled (virtual) oscillator and that of independent oscillators.

Another interesting thing we can see is that ${\mathfrak{R}}({b}_{{ii}}-({N}_{i}-1){b}_{i})$ has to be less than 0 to prevent the wave function from exploding at the boundary. Actually we haven't put any arbitrary constraint on the wave function. If ${\mathfrak{R}}({b}_{{ii}}-({N}_{i}-1){b}_{i})\geqslant 0$, even though it seems to be a 'good' wave function in the original coordinate, indeed the state will fail to satisfy the normalization condition if we really test it by integration.

With (25), we can then diagonalize the coupled X_i, finding their normal coordinates. Being independent from each other, we can then apply the propagator [55],

$$\begin{eqnarray}&&K({x}_{f},{x}_{i},t)=\sqrt{\displaystyle \frac{m\omega }{2\pi i{\rm{\hslash }}\mathrm{sin}(\omega t)}}\mathrm{exp}\left\{\displaystyle \frac{{im}\omega }{2{\rm{\hslash }}\mathrm{sin}(\omega t)}[({x}_{f}^{2}+{x}_{i}^{2})\mathrm{cos}(\omega t)-2{x}_{f}{x}_{i}]\right\}\end{eqnarray} \tag{ 32 }$$

on (31), so that

$$\begin{eqnarray}&&\psi (x,t)={\int }_{-\infty }^{\infty }K(x,x^{\prime} ,t)\psi (x^{\prime} ,0){\rm{}}\;{\rm{d}}x^{\prime} .\end{eqnarray} \tag{ 33 }$$

Transforming back to the ${\boldsymbol{X}}$ and ${\boldsymbol{u}}$ coordinates, now we've known how this system evolves.

Because here we have groups of identical particles, instead of really calculating the expectation value of a particular particle, there are many properties to exploit to bypass many procedures. In this system, the primary observables are position and momentum, so we will only consider them. As a reminder, ${{\boldsymbol{u}}}_{i}={{\rm{O}}}_{i}{{\boldsymbol{x}}}_{i}={({X}_{i},{u}_{i,2},{u}_{i,3},\ldots ,{u}_{i,{N}_{i}})}^{{\rm{T}}}$, and ${{\boldsymbol{P}}}_{i}={({P}_{i,1},{P}_{i,2},\ldots ,{P}_{i,{N}_{i}})}^{{\rm{T}}}={{\rm{O}}}_{i}{{\boldsymbol{p}}}_{i}={{\rm{O}}}_{i}{({p}_{i,1},{p}_{i,2},\ldots ,{p}_{i,{N}_{i}})}^{{\rm{T}}},$ where ${{\rm{O}}}_{i}$ is an orthogonal matrix. The particles in the original ${{\boldsymbol{x}}}_{i}$ coordinate are identical.

B.1.1. Single position operator

If we want to calculate the expectation value of ${x}_{i,j},$ since particles of this group are identical, we immediately recognize

$$\begin{eqnarray}&&\langle {x}_{i,j}\rangle =\langle {x}_{i,1}\rangle \;=\;\cdots \;=\;\langle {x}_{i,{N}_{i}}\rangle =\displaystyle \frac{\sqrt{{N}_{i}}}{{N}_{i}}\langle {X}_{i}\rangle =\displaystyle \frac{1}{\sqrt{{N}_{i}}}\langle {X}_{i}\rangle .\end{eqnarray} \tag{ 34 }$$

B.1.2. Single momentum operator

The momentum operator with respect to X_i is

$$\begin{eqnarray}&&{P}_{i,1}={[{({{\rm{O}}}_{i}^{-1})}^{{\rm{T}}}{{\boldsymbol{p}}}_{i}]}_{1}={({{\rm{O}}}_{i}{{\boldsymbol{p}}}_{i})}_{1}=\displaystyle \frac{1}{\sqrt{{N}_{i}}}\displaystyle \sum _{j=1}^{{N}_{i}}{p}_{i,j},\end{eqnarray} \tag{ 35 }$$

so we have

$$\begin{eqnarray}&&\langle {p}_{i,j}\rangle =\langle {p}_{i,1}\rangle \;=\;\cdots \;=\;\langle {p}_{i,{N}_{i}}\rangle =\displaystyle \frac{\sqrt{{N}_{i}}}{{N}_{i}}\langle {P}_{i,1}\rangle =\displaystyle \frac{1}{\sqrt{{N}_{i}}}\langle {P}_{i,1}\rangle .\end{eqnarray} \tag{ 36 }$$

B.1.3. Product of operators of the same group

Because we performed an orthogonal transformation,

$$\begin{eqnarray}&&\displaystyle \sum _{j=1}^{{N}_{i}}{x}_{i,j}^{2}={X}_{i}^{2}+\displaystyle \sum _{j=2}^{{N}_{i}}{u}_{i,j}^{2}\Rightarrow \langle {x}_{i,j}^{2}\rangle =\displaystyle \frac{1}{{N}_{i}}\left\langle {X}_{i}^{2}+\displaystyle \sum _{\alpha =2}^{{N}_{i}}{u}_{i,\alpha }^{2}\right\rangle \end{eqnarray} \tag{ 37 }$$

and

$$\begin{eqnarray}&&\langle {x}_{i,k}{x}_{i,l}\rangle =\displaystyle \frac{1}{{N}_{i}^{2}-{N}_{i}}\left\langle ({N}_{i}-1){X}_{i}^{2}-\displaystyle \sum _{\alpha =2}^{{N}_{i}}{u}_{i,\alpha }^{2}\right\rangle ,\ k\ne l.\end{eqnarray} \tag{ 38 }$$

Similarly we have

$$\begin{eqnarray}&&\displaystyle \sum _{j=1}^{{N}_{i}}{p}_{i,j}^{2}=\displaystyle \sum _{j=1}^{{N}_{i}}{P}_{i,j}^{2}\Rightarrow \langle {p}_{i,j}^{2}\rangle =\displaystyle \frac{1}{{N}_{i}}\left\langle \displaystyle \sum _{\alpha =1}^{{N}_{i}}{P}_{i,\alpha }^{2}\right\rangle ,\end{eqnarray} \tag{ 39 }$$

and

$$\begin{eqnarray}&&\langle {p}_{i,k}{p}_{i,l}\rangle =\displaystyle \frac{1}{{N}_{i}^{2}-{N}_{i}}\left\langle ({N}_{i}-1){P}_{i,1}^{2}-\displaystyle \sum _{\alpha =2}^{{N}_{i}}{P}_{i,\alpha }^{2}\right\rangle ,\ k\ne l.\end{eqnarray} \tag{ 40 }$$

A trickier one:

$$\begin{eqnarray} & & {{\boldsymbol{P}}}_{i}\cdot {{\boldsymbol{u}}}_{i}={{\boldsymbol{P}}}_{i}^{{\rm{T}}}{{\boldsymbol{u}}}_{i}={{\boldsymbol{p}}}_{i}^{{\rm{T}}}{{\rm{O}}}_{i}^{{\rm{T}}}{{\rm{O}}}_{i}{{\boldsymbol{x}}}_{i}={{\boldsymbol{p}}}_{i}\cdot {{\boldsymbol{x}}}_{i}\\ & \Rightarrow & \left\langle \displaystyle \sum _{j=1}^{{N}_{i}}{p}_{i,j}{x}_{i,j}\right\rangle =\left\langle \displaystyle \sum _{j=1}^{{N}_{i}}{P}_{i,j}{u}_{i,j}\right\rangle ,\end{eqnarray} \tag{ 41 }$$

with which we get

$$\begin{eqnarray}&&\langle {p}_{i,j}{x}_{i,j}\rangle =\displaystyle \frac{1}{{N}_{i}}\left\langle \left({X}_{i}{P}_{i,1}+\displaystyle \sum _{\alpha =2}^{{N}_{i}}{P}_{i,\alpha }{u}_{i,\alpha }\right)\right\rangle \end{eqnarray} \tag{ 42 }$$

and

$$\begin{eqnarray}&&\langle {p}_{i,j}{x}_{i,k}\rangle =\displaystyle \frac{1}{{N}_{i}^{2}-{N}_{i}}\left\langle ({N}_{i}-1){P}_{i,1}{X}_{i}-\displaystyle \sum _{\alpha =2}^{{N}_{i}}{P}_{i,\alpha }{u}_{i,\alpha }\right\rangle ,\ j\ne k.\end{eqnarray} \tag{ 43 }$$

B.1.4. Product of operators from different groups

Firstly,

$$\begin{eqnarray}&&\langle {x}_{i,k}{x}_{j,l}\rangle =\displaystyle \frac{1}{{N}_{i}{N}_{j}}\langle {X}_{i}{X}_{j}\rangle ,\ i\ne j,\end{eqnarray} \tag{ 44 }$$

and likely,

$$\begin{eqnarray}&&\left\{\begin{array}{l}\langle {p}_{i,k}{p}_{j,l}\rangle =\displaystyle \frac{1}{\sqrt{{N}_{i}{N}_{j}}}\langle {P}_{i,1}{P}_{j,1}\rangle \\ \langle {p}_{i,k}{x}_{j,l}\rangle =\displaystyle \frac{1}{\sqrt{{N}_{i}{N}_{j}}}\langle {P}_{i,1}{X}_{j}\rangle \end{array}\right.,\quad i\ne j.\end{eqnarray} \tag{ 45 }$$

If the coordinates ${u}_{i,j}$ for $j\gt 1$ are still identical after the transformation (21) like our models, then we can simplify our results derived above further. Below are the results. Note that these two matrices are symmetric and that the index j below is always larger than 1, and for this reason we replace ${u}_{i,j}$ by ${u}_{i}.$

The covariance matrix of two particles from different groups 1 and 2 is

$$\begin{eqnarray}\left(\begin{array}{cccc}\left\langle \tfrac{{X}_{1}^{2}}{{N}_{1}}+\tfrac{{N}_{1}-1}{{N}_{1}}{u}_{1}^{2}\right\rangle & {\mathfrak{R}}\left\langle \tfrac{{X}_{1}{P}_{\mathrm{1,1}}}{{N}_{1}}+\tfrac{{N}_{1}-1}{{N}_{1}}{u}_{1}{P}_{1,j}\right\rangle & {\mathfrak{R}}\left\langle \tfrac{{X}_{1}{X}_{2}}{\sqrt{{N}_{1}{N}_{2}}}\right\rangle & {\mathfrak{R}}\left\langle \tfrac{{X}_{1}{P}_{\mathrm{2,1}}}{\sqrt{{N}_{1}{N}_{2}}}\right\rangle \\ {\mathfrak{R}}\left\langle \tfrac{{X}_{1}{P}_{\mathrm{1,1}}}{{N}_{1}}+\tfrac{{N}_{1}-1}{{N}_{1}}{u}_{1}{P}_{1,j}\right\rangle & {\mathfrak{R}}\left\langle \tfrac{{P}_{1,1}^{2}}{{N}_{1}}+\tfrac{{N}_{1}-1}{{N}_{1}}{P}_{1,j}^{2}\right\rangle & {\mathfrak{R}}\left\langle \tfrac{{X}_{2}{P}_{\mathrm{1,1}}}{\sqrt{{N}_{1}{N}_{2}}}\right\rangle & {\mathfrak{R}}\left\langle \tfrac{{P}_{\mathrm{1,1}}{P}_{\mathrm{2,1}}}{\sqrt{{N}_{1}{N}_{2}}}\right\rangle \\ {\mathfrak{R}}\left\langle \tfrac{{X}_{1}{X}_{2}}{\sqrt{{N}_{1}{N}_{2}}}\right\rangle & {\mathfrak{R}}\left\langle \tfrac{{X}_{2}{P}_{\mathrm{1,1}}}{\sqrt{{N}_{1}{N}_{2}}}\right\rangle & \left\langle \tfrac{{X}_{2}^{2}}{{N}_{2}}+\tfrac{{N}_{2}-1}{{N}_{2}}{u}_{2}^{2}\right\rangle & {\mathfrak{R}}\left\langle \tfrac{{X}_{2}{P}_{\mathrm{2,1}}}{{N}_{2}}+\tfrac{{N}_{2}-1}{{N}_{2}}{u}_{2}{P}_{2,j}\right\rangle \\ {\mathfrak{R}}\left\langle \tfrac{{X}_{1}{P}_{\mathrm{2,1}}}{\sqrt{{N}_{1}{N}_{2}}}\right\rangle & {\mathfrak{R}}\left\langle \tfrac{{P}_{\mathrm{1,1}}{P}_{\mathrm{2,1}}}{\sqrt{{N}_{1}{N}_{2}}}\right\rangle & {\mathfrak{R}}\left\langle \tfrac{{X}_{2}{P}_{2}}{{N}_{2}}+\tfrac{{N}_{2}-1}{{N}_{2}}{u}_{2}{P}_{2,j}\right\rangle & \left\langle \tfrac{{P}_{2,1}^{2}}{{N}_{2}}+\tfrac{{N}_{2}-1}{{N}_{2}}{P}_{2,j}^{2}\right\rangle \end{array}\right)\end{eqnarray} \tag{ 46 }$$

and the covariance matrix of two particles both from the same group is

$$\begin{eqnarray}\left(\begin{array}{cccc}\left\langle \tfrac{{X}^{2}+(N-1){u}^{2}}{N}\right\rangle & {\mathfrak{R}}\left\langle \tfrac{{{XP}}_{1}+(N-1){{uP}}_{j}}{N}\right\rangle & \left\langle \tfrac{{X}^{2}-{u}^{2}}{N}\right\rangle & {\mathfrak{R}}\left\langle \tfrac{{{XP}}_{1}-{{uP}}_{j}}{N}\right\rangle \\ {\mathfrak{R}}\left\langle \tfrac{{{XP}}_{1}+(N-1){{uP}}_{j}}{N}\right\rangle & \left\langle \tfrac{{P}_{1}^{2}+(N-1){P}_{j}^{2}}{N}\right\rangle & {\mathfrak{R}}\left\langle \tfrac{{{XP}}_{1}-{{uP}}_{j}}{N}\right\rangle & \left\langle \tfrac{{P}_{1}^{2}-{P}_{j}^{2}}{N}\right\rangle \\ \left\langle \tfrac{{X}^{2}-{u}^{2}}{N}\right\rangle & {\mathfrak{R}}\left\langle \tfrac{{{XP}}_{1}-{{uP}}_{j}}{N}\right\rangle & \left\langle \tfrac{{X}^{2}+(N-1){u}^{2}}{N}\right\rangle & {\mathfrak{R}}\left\langle \tfrac{{{XP}}_{1}+(N-1){{uP}}_{j}}{N}\right\rangle \\ {\mathfrak{R}}\left\langle \tfrac{{{XP}}_{1}-{{uP}}_{j}}{N}\right\rangle & \left\langle \tfrac{{P}_{1}^{2}-{P}_{j}^{2}}{N}\right\rangle & {\mathfrak{R}}\left\langle \tfrac{{{XP}}_{1}+(N-1){{uP}}_{j}}{N}\right\rangle & \left\langle \tfrac{{P}_{1}^{2}+(N-1){P}_{j}^{2}}{N}\right\rangle \end{array}\right),\end{eqnarray} \tag{ 47 }$$

where for the sack of simplicity we ignore the first index i.

The state of this system is of this form

$$\begin{eqnarray}&&\psi ={ \mathcal N }\mathrm{exp}\left[-a\displaystyle \sum _{i=1}^{N}{x}_{i}^{2}+2b\displaystyle \sum _{i,j\gt i}^{N}{x}_{i}{x}_{j}\right],\end{eqnarray} \tag{ 48 }$$

a Gaussian state with permutation symmetry. The logarithmic negativity of Gaussian states is calculated with the methods shown in [38, 40], but in general the final result is too complex to give useful insights. However, if we assume that a and b in (48) are real, then it's possible to reduce it to such a simple form

$$\begin{eqnarray}{\mathsf{L}}=\left\{\begin{array}{ll}\tfrac{1}{2}\mathrm{ln}d & {\rm{if}}\;{\rm{}}d\gt 1\\ 0 & {\rm{otherswise}}\end{array}\right.,{\rm{}}\;{\rm{where}}\;{\rm{}}d=\left\{\begin{array}{ll}\tfrac{a+b}{a-b} & {\rm{if}}\;{\rm{}}b\geqslant 0\\ \tfrac{a+b-{bN}}{a+3b-{bN}} & {\rm{otherwise}}\end{array}\right.,\end{eqnarray} \tag{ 49 }$$

with basic but lengthy algebras, so we aren't going to show the procedure. The above equation clearly shows that when b = 0 the entanglement measure vanishes, which corresponds to the fact that at b = 0 the state is separable. Hence, we're more interested in circumstances where $b\ne 0$, so defining $r\equiv a/b$ we have

$$\begin{eqnarray}d=\left\{\begin{array}{ll}1+\tfrac{2}{r-1} & {\rm{if}}\;{\rm{}}b\geqslant 0\\ 1+\tfrac{2}{N-r-3} & {\rm{otherwise}}\end{array}\right..\end{eqnarray} \tag{ 50 }$$

There are also constraints for the wave function (48) such that it's square integrable:

$$\begin{eqnarray}a-(N-1)b\gt 0\;\& \;a+b\gt 0\;\& \;a\gt 0;{\rm{}}\;{\rm{or}}\;{\rm{}}\left\{\begin{array}{ll}r\gt N-1 & {\rm{if}}\;{\rm{}}b\geqslant 0\\ r\lt -1 & {\rm{otherwise}}\end{array}\right.,\end{eqnarray} \tag{ 51 }$$

taking into account $N\geqslant 2.$ Now let us discuss the value of d and thus ${\mathsf{L}}$ by the signs of b.

C.1.1.  $b\lt 0$

d in this case monotonically increases with r, except for the discontinuity at $N-r-3=0,$ namely at $r=N-3\geqslant -1$. Under the constraint (51), r never crosses this discontinuity, and thus the suprema of d and ${\mathsf{L}}$ are given by

$$\begin{eqnarray}&&\underset{r\to -{1}^{-}}{\mathrm{lim}}d=1+\displaystyle \frac{2}{N-2}{\rm{}}\;{\rm{and}}\;{\rm{}}\underset{r\to -{1}^{-}}{\mathrm{lim}}{\mathsf{L}}=\displaystyle \frac{1}{2}\mathrm{ln}\left(1+\displaystyle \frac{2}{N-2}\right),\end{eqnarray} \tag{ 52 }$$

because of the monotonicity of the logarithm. Now it's obvious that the maximum value of ${\mathsf{L}}$ decreases as N increases. In addition, at N = 2 both d and ${\mathsf{L}}$ aren't bounded from above.

At large N, because $r\lt -1,1/(N-r-3)\ll 1$ for all possible r and therefore we have

$$\begin{eqnarray}&&{\mathsf{L}}\simeq \displaystyle \frac{1}{2}\displaystyle \frac{2}{N-r-3}\simeq \displaystyle \frac{1}{N-r}.\end{eqnarray} \tag{ 53 }$$

Whenever $| r| \ll N,{\mathsf{L}}\simeq 1/N.$ Hence, the larger N is, the larger the state space where ${\mathsf{L}}\simeq 1/N$ is.

C.1.2.  $b\gt 0$

d now monotonically decreases as r increases. Since the discontinuity happens at r = 1, while the constraint (51) requires $r\gt N-1\geqslant 1$, the suprema of d and ${\mathsf{L}}$ are

$$\begin{eqnarray}&&\underset{r\to {(N-1)}^{+}}{\mathrm{lim}}d=1+\displaystyle \frac{2}{N-2}{\rm{}}\;{\rm{and}}\;{\rm{}}\underset{r\to {(N-1)}^{+}}{\mathrm{lim}}{\mathsf{L}}=\displaystyle \frac{1}{2}\mathrm{ln}\left(1+\displaystyle \frac{2}{N-2}\right),\end{eqnarray} \tag{ 54 }$$

same as the suprema for $b\lt 0,$ (52). And again at N = 2 there's no upper bound.

As (49) or (50) shows, d adopts different forms according to the sign of b. More interestingly while d for $b\gt 0$ doesn't depend on $N,d$ for $b\lt 0$ does. At first glance there seems to be discontinuity, but as (49) shows d under both circumstances approaches to 0 as b approaching 0, so the entanglement measure changes continuously as b switches the sign.

Now the state is

$$\begin{eqnarray}&&\psi ({{\boldsymbol{x}}}_{1},{{\boldsymbol{x}}}_{2})={ \mathcal N }\mathrm{exp}\left[-\displaystyle \sum _{i=1}^{{N}_{1}}{b}_{11}{x}_{1,i}^{2}+2{b}_{1}\displaystyle \sum _{i\gt j}^{{N}_{1}}{x}_{1,i}{x}_{1,j}-\displaystyle \sum _{i=1}^{{N}_{2}}{b}_{22}{x}_{2,i}^{2}+2{b}_{2}\displaystyle \sum _{i\gt j}^{{N}_{2}}{x}_{2,i}{x}_{2,j}+{c}_{12}{X}_{1}{X}_{2}\right],\end{eqnarray} \tag{ 55 }$$

where X₁ and X₂ are the normalized center of mass coordinates as defined in appendix A. We want to know the upper bound of ${{\mathsf{L}}}_{{ii}}$ of it. Because the two groups are of the same nature, the upper bound of them should have the same form, and we will just calculate ${{\mathsf{L}}}_{11}.$ Often it's easier to work in the normal coordinate,

$$\begin{eqnarray}&&\psi ={ \mathcal N }\mathrm{exp}[-{\epsilon }_{1}{X}_{1}^{2}-{\epsilon }_{2}{X}_{2}^{2}+{c}_{12}{X}_{1}{X}_{2}]\mathrm{exp}\left[-{\zeta }_{1}\displaystyle \sum _{i=1}^{{N}_{1}-1}{u}_{1,i}^{2}-{\zeta }_{2}\displaystyle \sum _{j=1}^{{N}_{2}-1}{u}_{2,j}^{2}\right],\end{eqnarray} \tag{ 56 }$$

where ${u}_{1,i}$ and ${u}_{2,j}$ are as defined in appendix A and ${\epsilon }_{i}\equiv {b}_{{ii}}-({N}_{i}-1){b}_{i}$ and ${\zeta }_{i}\equiv {b}_{{ii}}+{b}_{i}.$ Again assuming ${b}_{{ii}},{b}_{i},{c}_{12}$ are all real (so are ${\epsilon }_{i}$ and ${\zeta }_{i}$), the logarithmic negativity ${{\mathsf{L}}}_{11}$ of this state can be reduced to

$$\begin{eqnarray}{{\mathsf{L}}}_{11} & = & \left\{\begin{array}{ll}\displaystyle \frac{1}{2}\mathrm{ln}d & {\rm{if}}\;{\rm{}}d\gt 1\\ 0 & {\rm{otherswise}}\end{array}\right.,{\rm{}}\;{\rm{where}}\;{\rm{}}d=\left[1+\displaystyle \frac{1}{{N}_{1}}\left(-2+\displaystyle \frac{{\epsilon }_{1}}{{\zeta }_{1}}-\displaystyle \frac{4{\epsilon }_{2}{\zeta }_{1}}{{c}_{12}^{2}-4{\epsilon }_{1}{\epsilon }_{2}}\right)\right.\\ & & {\left.+\displaystyle \frac{\sqrt{{({\epsilon }_{1}{c}_{12}^{2}-4{\epsilon }_{1}^{2}{\epsilon }_{2}+4{\epsilon }_{2}{\zeta }_{1}^{2})}^{2}}}{{N}_{1}{\zeta }_{1}({c}_{12}^{2}-4{\epsilon }_{1}{\epsilon }_{2})}\right]}^{-1}\end{eqnarray} \tag{ 57 }$$

after basic but lengthy algebras. Therefore, d depends on the sign of

$$\begin{eqnarray}&&\kappa \equiv {\epsilon }_{1}{c}_{12}^{2}-4{\epsilon }_{1}^{2}{\epsilon }_{2}+4{\epsilon }_{2}{\zeta }_{1}^{2}.\end{eqnarray} \tag{ 58 }$$

C.2.1.  $\kappa \gt 0$

If it's larger than 0, then

$$\begin{eqnarray}&&d=\displaystyle \frac{{\zeta }_{1}{N}_{1}}{2{\epsilon }_{1}+({N}_{1}-2){\zeta }_{1}}=1+\displaystyle \frac{2{b}_{1}}{{b}_{11}-{b}_{1}}.\end{eqnarray} \tag{ 59 }$$

As before, there are constraint on these state parameters to make the wave function square integrable. Here, the constraints on b₁ and b₁₁ are

$$\begin{eqnarray}&&{b}_{11}\gt 0,{b}_{11}+{b}_{1}\gt 0{\rm{}}\;{\rm{and}}\;{\rm{}}{b}_{11}-({N}_{1}-1){b}_{1}\gt 0.\end{eqnarray} \tag{ 60 }$$

If ${b}_{1}\leqslant 0,$ then $d\leqslant 1$ because ${b}_{11}\gt 0$ and there's no entanglement. Since here we want to find the upper bound of ${{\mathsf{L}}}_{11},$ we don't have to take this into account. When ${b}_{1}\gt 0$, we only have to consider ${b}_{11}-{b}_{1}\gt 0$, because otherwise $d\lt 1$.

Now assuming ${b}_{11}-{b}_{1}\gt 0$ and ${b}_{1}\gt 0$ for the reasons stated above, the last constraint ${b}_{11}-({N}_{1}-1){b}_{1}\gt 0$ leads to

$$\begin{eqnarray}&&{b}_{11}-{b}_{1}\gt ({N}_{1}-2){b}_{1}\Rightarrow \displaystyle \frac{{b}_{1}}{({N}_{1}-2){b}_{1}}=\displaystyle \frac{1}{{N}_{1}-2}\gt \displaystyle \frac{{b}_{1}}{{b}_{11}-{b}_{1}},\end{eqnarray} \tag{ 61 }$$

and thus

$$\begin{eqnarray}&&d\lt 1+\displaystyle \frac{2}{{N}_{1}-2}\end{eqnarray} \tag{ 62 }$$

C.2.2.  $\kappa \lt 0$

If it's smaller than 0, then

$$\begin{eqnarray}&&d={\left(1-\displaystyle \frac{2}{{N}_{1}}+\displaystyle \frac{8{\epsilon }_{2}{\zeta }_{1}}{{N}_{1}(4{\epsilon }_{1}{\epsilon }_{2}-{c}_{12}^{2})}\right)}^{-1}.\end{eqnarray} \tag{ 63 }$$

To ensure the wave function is square integrable, we need the constraints

$$\begin{eqnarray}&&4{\epsilon }_{1}{\epsilon }_{2}-{c}_{12}^{2}\gt 0,{\epsilon }_{2}\gt 0{\rm{}}\;{\rm{and}}\;{\rm{}}{\zeta }_{1}\gt 0,\end{eqnarray} \tag{ 64 }$$

which implies

$$\begin{eqnarray}&&d\lt {\left(1-\displaystyle \frac{2}{{N}_{1}}\right)}^{-1}=1+\displaystyle \frac{2}{{N}_{1}-2}.\end{eqnarray} \tag{ 65 }$$

Because d has the same upper bound whether $\kappa \gt 0$ or $\kappa \lt 0$, we get

$$\begin{eqnarray}&&d\lt 1+\displaystyle \frac{2}{{N}_{1}-2}{\rm{}}\;{\rm{or}}\;{\rm{}}\mathrm{sup}{{\mathsf{L}}}_{11}=\displaystyle \frac{1}{2}\mathrm{ln}\left(1+\displaystyle \frac{2}{{N}_{1}-2}\right),\end{eqnarray} \tag{ 66 }$$

which has the same form as the one-group $\mathrm{sup}{\mathsf{L}}$ except N is replaced by ${N}_{1}.$

Let us begin with

$$\begin{eqnarray}&&\mathrm{sup}{\mathsf{L}}=\displaystyle \frac{1}{2}\mathrm{ln}\left(1+\displaystyle \frac{2}{N-2}\right)\simeq \displaystyle \frac{1}{N-2}\simeq \displaystyle \frac{1}{N}.\end{eqnarray} \tag{ 67 }$$

$1/N$ usually is a good approximation for $\mathrm{sup}{\mathsf{L}}$ for large N. But if we want an approximation that also serves as an upper bound, then $1/(N-2)$ should be chosen. The reason is

$$\begin{eqnarray}&&\displaystyle \frac{1}{N}\lt \mathrm{sup}{\mathsf{L}}\lt \displaystyle \frac{1}{N-2}.\end{eqnarray} \tag{ 68 }$$

The proof of the above inequality is simple. First,

$$\begin{eqnarray}&&\displaystyle \frac{d(1/N)}{{dN}}=-\displaystyle \frac{1}{{N}^{2}},\displaystyle \frac{d(\mathrm{sup}{\mathsf{L}})}{{dN}}=-\displaystyle \frac{1}{N(N-2)},\displaystyle \frac{d}{{dN}}\left(\displaystyle \frac{1}{N-2}\right)=-\displaystyle \frac{1}{{(N-2)}^{2}},\end{eqnarray} \tag{ 69 }$$

and thus

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dN}}\left(\displaystyle \frac{1}{N}-\mathrm{sup}{\mathsf{L}}\right)\gt 0{\rm{}}\;{\rm{and}}\;{\rm{}}\displaystyle \frac{d}{{dN}}\left(\mathrm{sup}{\mathsf{L}}-\displaystyle \frac{1}{N-2}\right)\gt 0{\rm{}}\;{\rm{for}}\;{\rm{}}N\gt 2,\end{eqnarray} \tag{ 70 }$$

meaning that both $\tfrac{1}{N}-\mathrm{sup}{\mathsf{L}}$ and $\mathrm{sup}{\mathsf{L}}-\tfrac{1}{N-2}$ monotonically increase for $N\gt 2$. Because

$$\begin{eqnarray}&&\underset{N\to \infty }{\mathrm{lim}}\displaystyle \frac{1}{N}=\underset{N\to \infty }{\mathrm{lim}}\mathrm{sup}{\mathsf{L}}=\underset{N\to \infty }{\mathrm{lim}}\displaystyle \frac{1}{N-2}=0,\end{eqnarray} \tag{ 71 }$$

$\tfrac{1}{N}-\mathrm{sup}{\mathsf{L}}$ and $\mathrm{sup}{\mathsf{L}}-\tfrac{1}{N-2}$ should approaches 0 as $N\to \infty .$ By their monotonicity, we have

$$\begin{eqnarray}&&\displaystyle \frac{1}{N}-\mathrm{sup}{\mathsf{L}}\lt 0{\rm{}}\;{\rm{and}}\;{\rm{}}\mathrm{sup}{\mathsf{L}}-\displaystyle \frac{1}{N-2}\lt 0.\end{eqnarray} \tag{ 72 }$$

This completes the proof. Therefore, ${\mathsf{L}}$ may be larger than $1/N$ but will always be less than $1/(N-2).$