format: text

text:
# Maximal entanglement versus entropy for mixed quantum states

## I Introduction

Over the past decade, the physical characteristic of the entanglement of quantum-mechanical states, both pure and mixed, has been recognized as a central resource in various aspects of quantum information processing. Significant settings include quantum communication Schumacher96 , cryptography Ekert91 , teleportation Bennett\_Brassard\_Crepeau\_Jozsa\_Peres\_Wootters , and, to an extent that is not quite so clear, quantum computation DiVincenzo95 . Given the central status of entanglement, the task of quantifying the degree to which a state is entangled is important for quantum information processing and, correspondingly, several measures of it have been proposed. These include: entanglement of formation Bennett\_DiVincenzo\_Smolin\_Wootters ; Wootters98 , entanglement of distillation Bennett\_Brassard\_Popescu\_Schumacher\_Smolin\_Wootters96 , relative entropy of entanglement Vedral\_Plenio\_Jacobs\_Knight , negativity Zyczkowski\_Horodecki\_Sanpera\_Lewenstein98 ; Vidal\_Werner02 , and so on. It is worth remarking that even for the smallest Hilbert space capable of exhibiting entanglement, i.e., the two-qubit system (for which Wootters has determined the entanglement of formation Wootters98 ), there are aspects of entanglement that remain to be explored.

Among the family of mixed quantum mechanical states, special status should be accorded to those that, for a given value of the entropy FOOT:entropy , have the largest possible degree of entanglement Munro\_James\_White\_Kwiat01 . The reason for this is that such states can be regarded as mixed-state generalizations of the Bell states, the latter being known to be the maximally entangled 2-qubit pure states. The notion of maximally entangled mixed states was introduced by Ishizaka and Hiroshima Ishizaka\_Hiroshima00  in a closely related setting, i.e., that of 2-qubit mixed states whose entanglement is maximized at fixed eigenvalues of the density matrix (rather than at fixed entropy of the density matrix). Evidently, the entanglement of the maximally entangled mixed states of Ishizaka and Hiroshima cannot be increased by any global unitary transformation. For these states, it was shown by Verstraete et al. Verstraete\_Audenaert\_deMoor01  that the maximality property continues to hold if any of the following three measures of entanglement—entanglement of formation, negativity, and relative entropy of entanglement—is replaced by one of the other two.

The question of the ordering of entanglement measures was raised by Eisert and Plenio Eisert\_Plenio99  and investigated numerically by them and analytically by Verstraete et al. Verstraete\_Audenaert\_Dehaene\_deMoor01 . It was proved by Virmani and Plenio Virmani\_Plenio00  that all good asymptotic entanglement measures are either identical or fail to uniformly give consistent orderings of density matrices. This implies that the resulting maximally entangled mixed states (MEMS) may depend on the measures one uses to quantify entanglement. Moreover, in finding the form of MEMS, one needs to quantify the mixedness of a state, and there can also be ordering problems for mixedness. This implies that the MEMS may depend on the measures of mixedness as well.

This Paper is organized as follows. We begin, in Secs. II and III, by reviewing several measures of entanglement and mixedness. In the main part of the Paper, Sec. IV, we consider various entanglement-versus-mixedness planes, in which entanglement and mixedness are quantified in several ways. Our primary objective, then, is to determine the frontiers, i.e., the boundaries of the regions occupied by physically allowed states in these planes, and to identify the structure of these maximally entangled mixed states. In Sec. V we make some concluding remarks.

## II Entanglement Criteria and their measures

It is well known that there are a large number of entanglement measures $`E`$. For a state described by the density matrix $`\rho `$ a good entanglement measure must satisfy, at least, the following conditions VedralPlenio98 ; Horodecki300 :
C1. (a) $`E(\rho )0`$; (b) $`E(\rho )=0`$ if $`\rho `$ is not entangled REF:SepEnt ; and (c) $`E(\text{Bell states})=1`$.
C2. For any state $`\rho `$ and any local unitary transformation, i.e., a unitary transformation of the form $`U_AU_B`$, the entanglement remains unchanged.
C3. Local operations, classical communication and postselection cannot increase the expectation value of the entanglement.
C4. Entanglement is convex under discarding information: $`_ip_iE(\rho _i)E(_ip_i\rho _i)`$.

The entanglement quantities chosen by us satisfy the properties C1–C4. Here, we do not impose the condition that any good entanglement measure should reduce to the entropy of entanglement (to be defined in the following) for pure states.

### II.1 Entanglement of formation and entanglement cost

The first measure we shall consider is the entanglement of formation $`E_\mathrm{F}`$ Bennett\_DiVincenzo\_Smolin\_Wootters ; it quantifies the amount of entanglement necessary to create the entangled state. It is defined by

$`E_\mathrm{F}(\rho )\underset{\{p_i,\psi _i\}}{\mathrm{min}}{\displaystyle \underset{i}{}}p_iE(|\psi _i\psi _i|),`$ (1)
where the minimization is taken over those probabilities $`\{p_i\}`$ and pure states $`\{\psi _i\}`$ that, taken together, reproduce the density matrix $`\rho =_ip_i|\psi _i\psi _i|`$. Furthermore, the quantity $`E(|\psi _i\psi _i|)`$ (usually called the entropy of entanglement) measures the entanglement of the pure state $`|\psi _i`$ and is defined to be the von Neumann entropy of the reduced density matrix $`\rho _i^{(A)}\text{Tr}_B|\psi _i\psi _i|`$, i.e.,

$`E(|\psi _i\psi _i|)=\text{Tr}\rho _i^{(A)}\mathrm{log}_2\rho _i^{(A)}.`$ (2)
For two-qubit systems, $`E_\mathrm{F}`$ can be expressed explicitly as Wootters98

$`E_\mathrm{F}(\rho )`$ $`=`$ $`h\left({\displaystyle \frac{1}{2}}[1+\sqrt{1C(\rho )^2}]\right),`$ (3a)
$`h(x)`$ $``$ $`x\mathrm{log}_2x(1x)\mathrm{log}_2(1x),`$ (3b)
where $`C(\rho )`$, the concurrence of the state $`\rho `$, is defined as
$`C(\rho )\mathrm{max}\{0,\sqrt{\lambda _1}\sqrt{\lambda _2}\sqrt{\lambda _3}\sqrt{\lambda _4}\},`$ (3c)
in which $`\lambda _1,\mathrm{},\lambda _4`$ are the eigenvalues of the matrix $`\rho (\sigma _y\sigma _y)\rho ^{}(\sigma _y\sigma _y)`$ in nondecreasing order and $`\sigma _y`$ is a Pauli spin matrix. $`E_\mathrm{F}(\rho )`$, $`C(\rho )`$, and the tangle $`\tau (\rho )C(\rho )^2`$ are equivalent measures of entanglement, inasmuch as they are monotonic functions of one another.

A measure associated with the entanglement of formation is the entanglement cost $`E_\mathrm{C}`$ Bennett\_DiVincenzo\_Smolin\_Wootters  which is defined via

$`E_\mathrm{C}(\rho )\underset{n\mathrm{}}{lim}{\displaystyle \frac{E_\mathrm{F}(\rho ^n)}{n}}.`$ (4)
This is the asymptotic value of the average entanglement of formation. $`E_\mathrm{C}`$ is, in general, difficult to calculate.

### II.2 Entanglement of distillation and relative entropy of entanglement

Related to the entanglement of formation is the entanglement of distillation $`E_\mathrm{D}`$ Bennett\_Brassard\_Popescu\_Schumacher\_Smolin\_Wootters96 , which characterizes the amount of entanglement of a state $`\rho `$ as the fraction of Bell states that can be distilled using the optimal purification procedure: $`E_\mathrm{D}(\rho )lim_n\mathrm{}m/n`$, where $`n`$ is the number of copies of $`\rho `$ used and $`m`$ is the maximal number of Bell states that can be distilled from them. The difference $`E_\mathrm{F}E_\mathrm{D}`$ can be regarded as undistillable entanglement. $`E_\mathrm{D}`$ is a difficult quantity to calculate, but the relative entropy of entanglement $`E_\mathrm{R}`$ Vedral\_Plenio\_Jacobs\_Knight , which we shall define shortly, provides an upper bound on $`E_\mathrm{D}`$ and is more readily calculable than it. For this reason, it is the second measure that we consider in this Paper. It is defined variationally via

$$E_\mathrm{R}(\rho )\underset{\sigma D}{\mathrm{min}}\mathrm{Tr}\left(\rho \mathrm{log}\rho \rho \mathrm{log}\sigma \right),$$
(5)
where $`D`$ represents the (convex) set of all separable density operators $`\sigma `$. In certain ways, the relative entropy of entanglement can be viewed as a distance $`𝒟(\rho ||\sigma ^{})`$ from the entangled state $`\rho `$ to the closest separable state $`\sigma ^{}`$. We remark that for pure states of two-qubit systems the relative entropy has the same value as the entanglement of formation.

### II.3 Negativity

The third measure that we shall consider is the negativity. The concept of the negativity of a state is closely related to the well-known Peres-Horodecki condition for the separability of a state Peres96\_Horodecki396 . If a state is separable (i.e., not entangled) then the partial transpose of its density matrix is again a valid state, i.e., it is positive semi-definite. It turns out that the partial transpose of a non-separable state has one or more negative eigenvalues. The negativity of a state Zyczkowski\_Horodecki\_Sanpera\_Lewenstein98  indicates the extent to which a state violates the positive partial transpose separability criterion. We will adopt the definition of negativity as twice the absolute value of the sum of the negative eigenvalues:

$$N(\rho )=2\mathrm{max}(0,\lambda _{\mathrm{neg}}),$$
(6)
where $`\lambda _{\mathrm{neg}}`$ is the sum of the negative eigenvalues of $`\rho ^{\mathrm{T}_B}`$. In $`C^2C^2`$ (i.e., two-qbit) systems it can be shown that the partial transpose of the density matrix can have at most one negative eigenvalue (see App. A). It was proved by Vidal and Werner Vidal\_Werner02  that negativity is an entanglement monotone, i.e., it satisfies criteria C1–C4 and, hence, is a good entanglement measure. We remark that for two-qubit pure states the negativity gives the same value as the concurrence does.

### II.4 Bures metric

The Bures metric of entanglement is defined as

$`E_\mathrm{B}(\rho )\underset{\sigma D}{\mathrm{min}}\left(22\sqrt{F(\rho ,\sigma )}\right),`$ (7)
where $`F(\rho ,\sigma )\left(\mathrm{Tr}\sqrt{\sqrt{\sigma }\rho \sqrt{\sigma }}\right)^2`$ is the fidelity. In the same way that relative entropy can, this entanglement measure can be viewed as the distance from the closest separable state to the entangled state considered, where the distance is now defined by $`𝒟(\rho ||\sigma )(22\sqrt{F(\rho ,\sigma )})`$ VedralPlenio98 . We remark that for two-qubit pure states the Bures metric reduces to the tangle defined in Sec. II.1.

### II.5 Lewenstein-Sanpera entanglement

It was shown by Lewenstein and Sanpera Lewenstein\_Sanpera98  that any density matrix $`\rho `$ has a decomposition into two parts:

$`\rho =\lambda \rho _\mathrm{s}+(1\lambda )\rho _\mathrm{e},`$ (8)
where $`\rho _\mathrm{s}`$ is separable, $`\rho _\mathrm{e}`$ is entangled, and the weight $`\lambda `$ is maximal, in which case the decomposition is unique. They refer to $`\rho _\mathrm{s}`$ as the best separable approximation (BSA) to $`\rho `$. It should be pointed out that, in general, it is nontrivial to establish the decomposition, even in the simplest relevant setting of $`C^2C^2`$ systems. Evidently, $`\lambda `$ and $`\rho _\mathrm{e}`$ contain information about the entanglement of $`\rho `$. Karnas and Lewenstein Karnas\_Lewenstein01  later showed that the quantity $`E_{\mathrm{LS}}(1\lambda )`$, which we will call the LS entanglement, satifies the above criteria, and hence is a good entanglement measure. For $`C^2C^2`$ systems it turns out that $`\rho _e`$ is a pure state, i.e., $`|\psi _\mathrm{e}\psi _\mathrm{e}|`$, and in this case it suggests that the quantity $`\stackrel{~}{E}_{\mathrm{LS}}`$, defined via

$`\stackrel{~}{E}_{\mathrm{LS}}(1\lambda )E(|\psi _\mathrm{e}|\psi _\mathrm{e}|),`$ (9)
may also be a good entanglement measure. We remark that the for two-qubit case the entanglement measure ($`1\lambda `$) is known to be equal to the Schmidt measure introduced in Ref. Eisert\_Briegel01 .

Even though the LS decomposition is not, in general, straightforward to find for the states in Eq. (22), the LS decomposition reads

$$\left(\begin{array}{cccc}x+\frac{r}{2}& 0& 0& \frac{r}{2}\\ 0& a& 0& 0\\ 0& 0& b& 0\\ \frac{r}{2}& 0& 0& y+\frac{r}{2}\end{array}\right)=\left(\begin{array}{cccc}x+\sqrt{ab}& 0& 0& \sqrt{ab}\\ 0& a& 0& 0\\ 0& 0& b& 0\\ \sqrt{ab}& 0& 0& y+\sqrt{ab}\end{array}\right)+(1\lambda )\rho _e,$$
(10)
where $`\rho _\mathrm{e}=|\varphi ^+\varphi ^+|`$ (with $`|\varphi ^\pm (|00\pm |11)/\sqrt{2}`$), $`(1\lambda )=r2\sqrt{ab}`$, and we consider only $`r2\sqrt{ab}>0`$ (as for $`r2\sqrt{ab}0`$ the whole density matrix is separable).

If we compute the concurrence of the state in Eq. (22), we find $`C=\mathrm{max}\{0,r2\sqrt{ab}\}`$. It is interesting to note that the LS entanglement gives, for the Ansatz states in Eq. (22), the same value as the concurrence does; but for an arbitrary state this is not, in general, true WellensKus01 .

### II.6 Ordering difficulties with entanglement measures

We now pause to touch on certain difficulties posed by the task of ordering physical states using entanglement. As first discussed and explored numerically by Eisert and Plenio Eisert\_Plenio99 , and subsequently investigated analytically by Verstraete et al. Verstraete\_Audenaert\_Dehaene\_deMoor01 , different entanglement measures can give different orderings for pairs of mixed states. This can be seen, e.g., from the plot of concurrence versus negativity, Fig. 1. The upper boundary is readily seen to be $`NC`$ whereas the lower boundary can be derived, giving $`N\sqrt{2(C\frac{1}{2})^2+\frac{1}{2}}+(C1)`$; see Ref. Verstraete\_Audenaert\_Dehaene\_deMoor01 .

Hence, when we wish to explain maximally entangled mixed states we need to be very explicit about the measure of entanglement (and also mixedness; see the following section). Different measures are likely to lead to different classes of MEMS states.

We end this section by mentioning the three entanglement measures that we shall use to compute the entanglement-mixedness frontiers: entanglement of formation, negativity, and relative entropy of entanglement. The first two of these are straightforward to compute, at least in two-qubit settings. For the third, certain results are available VedralPlenio98 ; Verstraete\_Audenaert\_deMoor01  that ease the computation of MEMS. We have also reviewed four additional measures (entanglement cost, entanglement of distillation, LS entanglement, and the Bures metric). Of these, however, the first two are rather difficult to compute, let alone maximize; the third is also difficult to compute, at least in practice. As for the fourth, calculating the entanglement involves finding the closest separable states (as is required for the case of relative entropy).

## III Measures of mixedness

In the entanglement-measure literature, two measures of mixedness have basically been used: $`\left(1\mathrm{Tr}[\rho ^2]\right)`$ and the von Neumann entropy. Whereas the latter has an natural significance stemming from its connections with statistical physics and information theory, the former is substantially easier to calculate. Of course, for density matrices that are almost completely mixed, the two measures show the same trend.

### III.1 von Neumann entropy

The von Neumann entropy, the standard measure of randomness of a statistical ensemble described by a density matrix, is defined by

$`S_\mathrm{V}(\rho )\mathrm{Tr}(\rho \mathrm{log}\rho )={\displaystyle \underset{i}{}}\lambda _i\mathrm{log}\lambda _i,`$ (11)
where $`\lambda _i`$ are the eigenvalues of the density matrix $`\rho `$ and the $`\mathrm{log}`$ is taken to base $`𝒩`$, the dimension of the Hilbert space in question. It is straightforward to show that the extremal values of $`S_\mathrm{V}`$ are zero (for pure states) and unity (for completely mixed states). To compute the von Neumann entropy it is necessary to have the full knowledge of the eigenvalue spectrum.

As we shall mention in the following subsection, there is a linear entropy threshold above which all states are separable. Qualitatively identical behavior is encountered for the von Neumann entropy. In particular, as we shall see in Sec. IV.3.1, for two-qubit systems all states are separable for $`S_\mathrm{V}(1/2)\mathrm{log}_4(1/12)0.896`$.

### III.2 Purity and linear entropy

The second measure that we shall consider is called the linear entropy and is based on the purity of a state, $`𝒫\mathrm{Tr}\left[\rho ^2\right]`$, which ranges from 1 (for a pure state) to $`1/𝒩`$ for a completely mixed state with dimension $`𝒩`$. The linear entropy $`S_\mathrm{L}`$ is defined via

$`S_\mathrm{L}(\rho ){\displaystyle \frac{𝒩}{𝒩1}}(1\text{Tr}[\rho ^2]),`$ (12)
which ranges from 0 (for a pure state) to 1 (for a maximally mixed state). The linear entropy is generally a simpler quantity to calculate than the von Neumann entropy as there is no need for diagonalization. For $`C^2C^2`$ systems the linear entropy can be written explicity as

$`S_\mathrm{L}(\rho ){\displaystyle \frac{4}{3}}(1\text{Tr}\left[\rho ^2\right]).`$ (13)
A related measure, which we shall not use in this Paper (but mention for the sake of completeness), is the inverse participation ratio. Defined via $`R1/\text{Tr}\left[\rho ^2\right]`$, it ranges from $`1`$ (for a pure state) to $`𝒩`$ (for the maximally mixed state). An attractive property of the inverse participation ratio is that all states with $`R𝒩1`$ are separable Zyczkowski\_Horodecki\_Sanpera\_Lewenstein98 , which implies all states with a linear entropy $`S_\mathrm{L}(\rho )𝒩(𝒩2)/\left(𝒩1\right)^2`$ (which is $`8/9`$ when $`𝒩=4`$) are separable.

### III.3 Comparing linear and von Neumann entropies

The aim of this subsection is to illustrate the difference between the linear and von Neumann entropies. We shall do this by considering the $`𝒩=4`$ Hilbert space, and seeking the highest and lowest von Neumann entropies consistent with a given value of linear entropy. Before restricting $`𝒩`$ to 4, the corresponding stationarity problem reads:

$`\delta \left(S_\mathrm{V}(\rho )+\beta {\displaystyle \frac{𝒩1}{2𝒩}}S_\mathrm{L}(\rho )(\nu 1)\text{Tr}\rho \right)=0,`$ (14)
where $`\beta `$ and $`\nu `$ are, respectively, Langrange multipliers that enforce the constraints that linear entropy be fixed and that $`\rho `$ be normalized. Thus, we arrive at the engaging self-consistency condition

$`\rho =\mathrm{exp}(\nu \beta \rho ),`$ (15)
in which $`\nu `$ and $`\beta `$ can be fixed upon implementing the constraints.

By hypothesizing certain relationships among the eigenvalues of the density matrix, we have been able, for the case of $`𝒩=4`$, to find certain solutions of the self-consistency condition and, hence, a candidate for the boundary of the region of $`S_\mathrm{L}`$ vs. $`S_\mathrm{V}`$ plane (see Fig. 2) corresponding to legitimate density matrices. These solutions, when given in terms of eigenvalues, read

$`\{{\displaystyle \frac{1r}{4}},{\displaystyle \frac{1r}{4}},{\displaystyle \frac{1r}{4}},{\displaystyle \frac{1+3r}{4}}\},\text{for}0r1,`$ (16a)
$`\{r,1r,0,0\},\text{for}{\displaystyle \frac{1}{2}}r1,`$ (16b)
$`\{r,{\displaystyle \frac{1r}{2}},{\displaystyle \frac{1r}{2}},0\},\text{for}0r{\displaystyle \frac{1}{3}},`$ (16c)
$`\{{\displaystyle \frac{4r}{12}},{\displaystyle \frac{4r}{12}},{\displaystyle \frac{4r}{12}},{\displaystyle \frac{3r}{12}}\},\text{for}0r1,`$ (16d)
and they correspond to the upper boundary, and the lowest, middle, and highest pieces of the lower boundary, respectively. Note that the lower boundary comprises three (in general, $`𝒩1`$) segments that meet at cusps. We remark, parenthetically, that the solutions with zero eigenvalues correspond to extrema within some subspace spanned by those eigenvectors with nonzero eigenvalues, and therefore only obey the stationarity condition (15) within the subspace.

Is there any significance to the boundary states? Boundary segment (a) includes the Werner states defined in Eq. (24). Boundary segment (b) includes the first branch of the MEMS for $`E_\mathrm{F}`$ and $`S_\mathrm{L}`$ specified below in Eq. (23). The segment (c) includes the states

$`\rho _c=r|\varphi ^+\varphi ^+|+{\displaystyle \frac{1r}{2}}\left(|0101|+|1010|\right).`$ (17)
States on segment (d) are all unentangled. Of course, the boundary segments include not only the specified states but also all states derivable from them by global unitary transformation.

As for the interior, we have obtained this numerically by constructing a large number of random sets of eigenvalues of legitimate density matrices, and computing for each the two entropies. As Fig. 2 shows, no points lie outside the boundary curve, providing confirmatory evidence for the forms given in Eq. (16).

The fact that the bounded region is two-dimensional indicates the lack of precision with which the linear entropy characterizes the von Neumann entropy (and vice versa, if one wishes). In particular, the figure reveals an ordering difficulty: pairs of states, $`A`$ and $`B`$, exist for which $`S_\mathrm{L}^AS_\mathrm{L}^B`$ and $`S_\mathrm{V}^AS_\mathrm{V}^B`$ differ in sign. Worse still, states having a common value of $`S_\mathrm{V}`$ have a continuum of values of $`S_\mathrm{L}`$, and vice versa.

## IV Entanglement-versus-mixedness frontiers

We now attempt to identify regions in the plane spanned by entanglement and mixedness that are inhabited by physical states (i.e., characterized by legitimate density matrices). We shall consider the various measures of entanglement and mixedness discussed in the previous section. Of particular interest will be the structure of the states that inhabit the frontier, i.e., the boundary delimiting the region of physical states. Frontier states are maximal in the following sense: for a given value of mixedness they are maximally entangled; for a given value of entanglement they are maximally mixed.

### IV.1 Parametrization of maximal states

The aim of this subsection is to derive the general form of the maximal states given in Eq. (21), which is what we will use to parametrize maximal states. In Ref. Verstraete\_Audenaert\_deMoor01 , it is shown that, given a fixed set of eigenvalues, all states that maximize one of the three entanglement measures (entanglement of formation, negativity or relative entropy) automatically maximize the other two. It was further shown that the global unitary transformation that takes arbitrary states into maximal ones has the form

$`U=(U_1U_2)TD_\varphi \mathrm{\Phi }^{},`$ (18)
where $`U_1`$ and $`U_2`$ are arbitary local unitary transformations,

$`T\left(\begin{array}{cccc}0& 0& 0& 1\\ \frac{1}{\sqrt{2}}& 0& \frac{1}{\sqrt{2}}& 0\\ \frac{1}{\sqrt{2}}& 0& \frac{1}{\sqrt{2}}& 0\\ 0& 1& 0& 0\end{array}\right),`$ (19)
$`D_\varphi `$ is a unitary diagonal matrix, and $`\mathrm{\Phi }`$ is the unitary matrix that diagonalizes the density matrix $`\rho `$, i.e., $`\rho =\mathrm{\Phi }\mathrm{\Lambda }\mathrm{\Phi }^{}`$, where $`\mathrm{\Lambda }`$ is a diagonal matrix, the diagonal elements of which are the four eigenvalues of $`\rho `$ listed in nonincreasing order ($`\lambda _1\lambda _2\lambda _3\lambda _4`$). Hence, the general form of a density matrix that is maximal, given a set of eigenvalues, is (up to local unitary transformations)

$$T\left(\begin{array}{cccc}\lambda _1& 0& 0& 0\\ 0& \lambda _2& 0& 0\\ 0& 0& \lambda _3& 0\\ 0& 0& 0& \lambda _4\end{array}\right)T^{}=\left(\begin{array}{cccc}\lambda _4& 0& 0& 0\\ 0& \frac{\lambda _1+\lambda _3}{2}& \frac{\lambda _1\lambda _3}{2}& 0\\ 0& \frac{\lambda _1\lambda _3}{2}& \frac{\lambda _1+\lambda _3}{2}& 0\\ 0& 0& 0& \lambda _2\end{array}\right).$$
(20)
This matrix is locally equivalent to the form

$$\left(\begin{array}{cccc}x+\frac{r}{2}& 0& 0& \frac{r}{2}\\ 0& a& 0& 0\\ 0& 0& b& 0\\ \frac{r}{2}& 0& 0& x+\frac{r}{2}\end{array}\right),$$
(21)
with $`x+\frac{r}{2}=(\lambda _1+\lambda _3)/2`$, $`r=\lambda _1\lambda _3`$, $`a=\lambda _2`$, and $`b=\lambda _4`$. The above derivation justifies the Ansatz form (22) used in Ref. Munro\_James\_White\_Kwiat01  to derive the entanglement of formation vs. linear entropy MEMS. We remark that one may as well use the four eigenvalues $`\lambda _i`$’s as the parametrization. Nevertheless, the form (21), as well as (22), can be nicely viewed as a mixture of a Bell state $`|\varphi ^+`$ with some diagonal separable mixed state.

### IV.2 Entanglement-versus-linear-entropy frontiers

We begin by measuring mixedness in terms of the linear entropy, and comparing the frontier states for various measures of entanglement.

#### IV.2.1 Entanglement of formation

The characterization of physical states in terms of their entanglement of formation and linear entropy was introduced by Munro et al. in Ref. Munro\_James\_White\_Kwiat01 . (Strictly speaking, they considered the tangle rather than the equivalent entanglement of formation.) Here, we shall consider yet another equivalent quantity: concurrence (see Sec. II.1). In order to find the frontier, Munro et al. proposed Ansatz states of the form

$`\rho _{\text{Ansatz}}=\left(\begin{array}{cccc}x+\frac{r}{2}& 0& 0& \frac{r}{2}\\ 0& a& 0& 0\\ 0& 0& b& 0\\ \frac{r}{2}& 0& 0& y+\frac{r}{2}\end{array}\right),`$ (22)
where $`x,y,a,b,r0`$ and $`x+y+a+b+r=1`$. They found that, of these, the subset

$`\rho _{\text{MEMS:}E_\mathrm{F},S_\mathrm{L}}=\{\begin{array}{cc}\rho _\mathrm{I}(r),\hfill & \mathrm{for}\frac{2}{3}r1;\hfill \\ \rho _{\mathrm{I}\mathrm{I}}(r),\hfill & \mathrm{for}0r\frac{2}{3};\hfill \end{array}`$ (23a)
$`\rho _\mathrm{I}(r)=\left(\begin{array}{cccc}\frac{r}{2}& 0& 0& \frac{r}{2}\\ 0& 1r& 0& 0\\ 0& 0& 0& 0\\ \frac{r}{2}& 0& 0& \frac{r}{2}\end{array}\right),\rho _{\mathrm{I}\mathrm{I}}(r)=\left(\begin{array}{cccc}\frac{1}{3}& 0& 0& \frac{r}{2}\\ 0& \frac{1}{3}& 0& 0\\ 0& 0& 0& 0\\ \frac{r}{2}& 0& 0& \frac{1}{3}\end{array}\right),`$ (23b)
lies on the boundary in the tangle vs. linear-entropy plane and, accordingly, named these MEMS, in the sense that these states have maximal tangle for a given linear entropy. We remark that at the crossing point of the two branches, $`r=2/3`$, the density matrices on either side coincide.

In Fig. 3 we plot the entanglement of formation/concurrence vs. linear entropy for the family of MEMS (23); this gives the frontier curve. For the sake of comparison, we also give the curve associated with the family of Werner states of the form

$$\rho _\mathrm{W}r|\varphi ^+\varphi ^+|+\frac{1r}{4}𝟙=\left(\begin{array}{cccc}\frac{\mathrm{𝟙}+𝕣}{\mathrm{𝟜}}& \mathrm{𝟘}& \mathrm{𝟘}& \frac{𝕣}{\mathrm{𝟚}}\\ \mathrm{𝟘}& \frac{\mathrm{𝟙}𝕣}{\mathrm{𝟜}}& \mathrm{𝟘}& \mathrm{𝟘}\\ \mathrm{𝟘}& \mathrm{𝟘}& \frac{\mathrm{𝟙}𝕣}{\mathrm{𝟜}}& \mathrm{𝟘}\\ \frac{𝕣}{\mathrm{𝟚}}& \mathrm{𝟘}& \mathrm{𝟘}& \frac{\mathrm{𝟙}+𝕣}{\mathrm{𝟜}}\end{array}\right).$$
(24)
Evidently, for a given value of linear entropy these MEMS (which we shall denote by $`\{\mathrm{MEMS}:E_\mathrm{F},S_\mathrm{L}\}`$) achieve the highest concurrence. As the tangle $`\tau `$ and entanglement of formation $`E_\mathrm{F}`$ are monotonic functions of the concurrence, Eq. (23) also gives the boundary curve for these measures. This raises an interesting question: Is (23) optimal for other measures of entanglement?

#### IV.2.2 Relative entropy as the entanglement measure

To find the frontier states for the relative entropy of entanglement we again turn our attention to the maximal density matrix (21). For this form of density matrix the linear entropy is given (with $`x`$ expressed in terms of $`a,b,r`$) by

$`S_\mathrm{L}={\displaystyle \frac{2}{3}}(3a^2+2a\left(1b\right)+\left(1b\right)\left(1+3b\right)r^2).`$ (25)
To calculate the relative entropy of entanglement, we need to determine the closest separable state to (21). It is simpler to do this analysis via several cases. We begin by considering the Rank-2 and Rank-3 cases of (21). We set $`b=0`$ ($`\lambda _4=0`$) and express $`x`$ in terms of $`a`$ and $`r`$ in the density matrix, obtaining

$`\rho =\left(\begin{array}{cccc}\frac{1a}{2}& 0& 0& \frac{r}{2}\\ 0& a& 0& 0\\ 0& 0& 0& 0\\ \frac{r}{2}& 0& 0& \frac{1a}{2}\end{array}\right),`$ (26)
and thus find that the closest separable density matrix $`\sigma ^{}`$ is given by VedralPlenio98

$`\sigma ^{}`$ $`=`$ $`\left(\begin{array}{cccc}C& 0& 0& D\\ 0& E& 0& 0\\ 0& 0& 12CE& 0\\ D& 0& 0& C\end{array}\right),`$ (27a)
$`C`$ $``$ $`{\displaystyle \frac{\left(1+a\right)\left(1a^2r^2\right)}{2\left(1+ar\right)\left(1+a+r\right)}},`$ (27b)
$`D`$ $``$ $`{\displaystyle \frac{a\left(1+a\right)r}{\left(1+ar\right)\left(1+a+r\right)}},`$ (27c)
$`E`$ $``$ $`{\displaystyle \frac{a\left(1+a\right)^2}{\left(1+ar\right)\left(1+a+r\right)}}.`$ (27d)
The relative entropy of entanglement is now simply given by

$`E_\mathrm{R}(\rho )=\mathrm{Tr}\left(\rho \mathrm{log}\rho \rho \mathrm{log}\sigma ^{}\right)`$
$`={\displaystyle \frac{1+a}{2}}\mathrm{log}{\displaystyle \frac{(1+a)^2r^2}{(1+a)^2}}+{\displaystyle \frac{r}{2}}\mathrm{log}{\displaystyle \frac{1+a+r}{1+ar}},`$ (28)
with the linear entropy being given by

$`S_\mathrm{L}={\displaystyle \frac{2}{3}}\left(1+2a3a^2r^2\right),`$ (29)
subject to the constraint $`(a+r)1`$. For the Rank-2 case $`a=1r`$ ($`b=x=0`$), and the resulting solution is the Rank-2 matrix $`\rho _\mathrm{I}(r)`$ given in Eq. (23) with $`1/2r1`$. We remark that this Rank-2 solution is always a candidate MEMS for the three entanglement measures that we consider in this Paper. In order to determine whether or in what range the Rank-2 solution achieves the global maximum, we need to compare it with the Rank-3 and Rank-4 solutions.

By maximizing $`E_\mathrm{R}(\rho )`$ for a given value of $`S_\mathrm{L}`$, we find the following stationary condition:

$`r\mathrm{log}{\displaystyle \frac{(1+a)^2r^2}{(1+a)^2}}=(3a1)\mathrm{log}{\displaystyle \frac{1+a+r}{1+ar}}.`$ (30)
Given a value of $`S_\mathrm{L}`$, we can solve Eqs. (29) and (30), at least numerically, to obtain the parameters $`a`$ and $`r`$, and hence, from Eq. (26), the Rank-3 MEMS. However, if the constraint inequality $`a+r1`$ turns out to be violated, the solution is invalid.

We now turn to the Rank-4 case. It is straightforward, if tedious, to show that the Werner states, Eq. (24), obey the stationarity conditions appropriate for Rank 4. However, it turns out that this solution is not maximal.

To summarize, the frontier states, which we denote by $`\{\mathrm{MEMS}:E_\mathrm{R},S_\mathrm{L}\}`$, are states of the form (26); the dependence of the parameters $`a`$ and $`r`$ on $`S_\mathrm{L}`$ is shown in Fig. 4. In Fig. 5, we show the resulting frontier, as well as curves corresponding to non-maximal stationary states. The frontier states have the following structure: (i) for $`S_\mathrm{L}0.5054`$ they are the Rank-2 MEMS of Eq. (23) but with $`r`$ restricted to the range from $`1`$ (at $`S_\mathrm{L}=0`$) to approximately $`0.7459`$ (at $`S_\mathrm{L}0.5054`$); (ii) for $`S_\mathrm{L}0.5054`$ the MEMS are Rank 3, with parameters $`a`$ and $`r`$ satisfying Eqs. (30) and (29) at each value of $`S_\mathrm{L}`$, and $`(a,r)`$ ranging between approximately $`(0.3056,0.7459)`$ (at $`S_\mathrm{L}0.5054`$) and $`(1/3,0)`$ (at $`S_\mathrm{L}=8/9`$). As noted previously, beyond $`S_\mathrm{L}=8/9`$ there are no entangled states. As the inset of Fig. 4 shows, the parameter $`a`$ can be regarded as a continous function of parameter $`r[0,1]`$. The two branches of the solution, (i) and (ii), cross at $`(S_\mathrm{L}^{},E_\mathrm{R}^{})(0.5054,0.3422)`$; at this point, the states on the two branches coincide,

$`\rho ^{}\left(\begin{array}{cccc}0.372947& 0& 0& 0.372947\\ 0& 0.254106& 0& 0\\ 0& 0& 0& 0\\ 0.372947& 0& 0& 0.372947\end{array}\right).`$ (31)
Just as in the case of entanglement of formation vs. linear entropy, the density matrix is continuous at the transition between branches.

We remark that the curve generated by the states $`\{\mathrm{MEMS}:E_\mathrm{F},S_\mathrm{L}\}`$, when plotted on the $`E_\mathrm{R}`$ vs. $`S_\mathrm{L}`$ plane, falls just slightly below that generated by the states $`\{\mathrm{MEMS}:E_\mathrm{R},S_\mathrm{L}\}`$ for $`S_\mathrm{L}0.5054`$ (and coincides for smaller values of $`S_\mathrm{L}`$).

We also remark that the parameter $`r`$ turns out to be the concurrence $`C`$ of the states, so that Fig. 4 can be interpreted a plot of the concurrence of the frontier states vs. their linear entropy. By comparing this concurrence vs. linear entropy curve to that in Fig. 3, we find that the former lies just slightly below the latter for $`S_\mathrm{L}0.5054`$ (and the two coincide for smaller values of $`S_\mathrm{L}`$), the maximal difference between the two being less than $`10^2`$.

It is evident that, for a given linear entropy, the relative entropies of entanglement for both $`\{\mathrm{MEMS}:E_\mathrm{R},S_\mathrm{L}\}`$ and $`\{\mathrm{MEMS}:E_\mathrm{F},S_\mathrm{L}\}`$ are significantly less than the corresponding entanglements of formation. In fact, for small degrees of impurity, the entanglements of formation for the two MEMS states are quite flat; however, the relative entropies of entanglement fall quite rapidly. More specifically, for a change in linear entropy of $`\mathrm{\Delta }S_\mathrm{L}=0.1`$ near $`S_\mathrm{L}=0`$ we have that $`\mathrm{\Delta }E_\mathrm{F}0.05`$ (see Fig. 3) and $`\mathrm{\Delta }E_R0.2`$ (see Fig. 5). As the curves of the states $`\{\mathrm{MEMS}:E_\mathrm{F},S_\mathrm{L}\}`$ and $`\{\mathrm{MEMS}:E_\mathrm{R},S_\mathrm{L}\}`$ are very close on the two planes, $`E_\mathrm{F}`$ vs. $`S_\mathrm{L}`$ and $`E_\mathrm{R}`$ vs. $`S_\mathrm{L}`$, we show, in Fig. 6, the entanglement difference $`E_\mathrm{F}E_\mathrm{R}`$ for the states $`\{\mathrm{MEMS}:E_\mathrm{F},S_\mathrm{L}\}`$, and compare it with the corresponding difference for the Werner states. While it is clear that $`E_\mathrm{R}(\rho )E_\mathrm{F}(\rho )`$, for certain values of the linear entropy the difference turns out to be quite large, this difference being uniformly larger for $`\{\mathrm{MEMS}:E_\mathrm{F},S_\mathrm{L}\}`$ than for the Werner state; see Fig. 6. This raises several interesting questions, e.g., whether, for a given value of $`E_\mathrm{F}`$, the undistillable entanglement (the difference between the entanglement of formation and the relative entropy of entanglement is the lower bound) increases as states become more mixed.

As we have seen, Werner states are not frontier states either in the case of entanglement of formation or in the case of relative entropy of entanglement. By contrast, as we shall see in the next section, if we measure entanglement via negativity, then for a given amount of linear entropy, the Werner states (as well as another Rank-3 class of states) achieve the largest value of entanglement. Said equivalently, the Werner states belong to $`\{\mathrm{MEMS}:\mathrm{N},\mathrm{S}_\mathrm{L}\}`$.

#### IV.2.3 Negativity

In order to derive the form of the MEMS in the case of negativity, we again consider the density matrix of the form (21), for which it is straightforward to show that the negativity $`N`$ is given by

$`N=\mathrm{max}\{0,\sqrt{(ab)^2+r^2}(a+b)\}.`$ (32)
Furthermore, because we aim to find the entanglement frontier, we can simply restrict our attention to states satisfying $`N>0`$, i.e., to states that are entangled Peres96\_Horodecki396 . Then, by making $`N`$ stationary at fixed $`S_\mathrm{L}`$ and with the constraint $`2x+a+b+r=1`$, we find two one-parameter families of stationary states (in addition to the Rank-2 MEMS, which are common to all three entanglement measures). The parameters of the first family obey

$`a=b=x,r=14x.`$ (33)
When expressed in terms of parameter $`r`$, the density matrix takes the form

$`\rho _{\text{MEMS:}N,S_\mathrm{L}}^{(1)}=\left(\begin{array}{cccc}\frac{1+r}{4}& 0& 0& \frac{r}{2}\\ 0& \frac{1r}{4}& 0& 0\\ 0& 0& \frac{1r}{4}& 0\\ \frac{r}{2}& 0& 0& \frac{1+r}{4}\end{array}\right),`$ (34)
which are precisely the Werner states in Eq. (24). For the second solution, the parameters obey

$$a=\frac{42\sqrt{3r^2+1}}{6},b=0,x=\frac{1+\sqrt{3r^2+1}}{6}\frac{r}{2}.$$
(35)
When expressed in terms of parameter $`r`$, the density matrix takes the form

$$\rho _{\text{MEMS:}N,S_\mathrm{L}}^{(2)}=\left(\begin{array}{cccc}\frac{1+\sqrt{3r^2+1}}{6}& 0& 0& \frac{r}{2}\\ 0& \frac{42\sqrt{3r^2+1}}{6}& 0& 0\\ 0& 0& 0& 0\\ \frac{r}{2}& 0& 0& \frac{1+\sqrt{3r^2+1}}{6}\end{array}\right).$$
(36)
We remark that the two solutions give the same bound on the negativity for a given value of linear entropy. The resulting frontier in the negativity vs. linear entropy plane is shown in Fig. 7.

We have shown by numerical exploration FOOT:NSL  that the family states found in the previous paragraph do indeed achieve the maximal negativity. The results of this exploration are indicated by points in the inset in Fig. 7, and compelling evidence for the maximality of the negativity is furnished by the fact that no points lie beyond the frontier.

Thus, the states $`\{\mathrm{MEMS}:N,S_\mathrm{L}\}`$ on the boundary include, up to local unitary transformations, both Werner states in Eq. (34) and states in Eq. (36). We also plot in Fig. 7 the curve belonging to $`\{\mathrm{MEMS}:E_\mathrm{F},S_\mathrm{L}\}`$; note that it falls slightly below the curve associated with $`\{\mathrm{MEMS}:N,S_\mathrm{L}\}`$ and that it has a cusp, due to the structure of the states, at the value $`2/3`$ for the parameter $`r`$ in Eq. (23). Here, we see that maximally entangled mixed states change their form when we adopt a different entanglement measure.

### IV.3 Entanglement-versus-von-Neumann-entropy frontiers

We continue this section by choosing to measure mixedness in terms of the von Neumann entropy, and comparing the frontier states for various measures of entanglement.

#### IV.3.1 Entanglement of formation

To find this frontier, we consider states of the form (21), and compute for them the concurrence and the von Neumann entropy:

$`C`$ $`=`$ $`r2\sqrt{ab},`$ (37a)
$`S_\mathrm{V}`$ $`=`$ $`a\mathrm{log}ab\mathrm{log}bx\mathrm{log}x(x+r)\mathrm{log}(x+r).`$ (37b)
Note that the parameters obey the normalization constraint $`2x+a+b+r=1`$.

As we remarked previously, the Rank 2 MEMS is always a candidate. For the Rank-3 case, we can set $`b=0`$ in Eq. (37). By maximizing $`C`$ at fixed $`S_\mathrm{V}`$, we find a stationary solution:
(i) $`r=C`$, $`x=(43C\sqrt{43C^2})/6`$, and $`a=(\sqrt{43C^2}1)/3`$; the resulting density matrix is

$`\rho _\mathrm{i}=\left(\begin{array}{cccc}\frac{4\sqrt{43C^2}}{6}& 0& 0& \frac{C}{2}\\ 0& \frac{\sqrt{43C^2}1}{3}& 0& 0\\ 0& 0& 0& 0\\ \frac{C}{2}& 0& 0& \frac{4\sqrt{43C^2}}{6}\end{array}\right).`$ (38)
For the Rank-4 case ($`b0`$), the stationarity condition can be shown to be

$`u\mathrm{log}(u)`$ $`=`$ $`w\mathrm{log}(w),`$ (39a)
$`2u\mathrm{log}(u)`$ $`=`$ $`(u+w)\mathrm{log}(v),`$ (39b)
where $`u\sqrt{a/(x+r)}`$, $`v\sqrt{x/(x+r)}`$, and $`w\sqrt{b/(x+r)}`$. There are two solutions, due to the two-to-one property of the function $`z\mathrm{log}z`$ for $`z(0,1)`$. The first one is ($`u=v=w`$).
(ii) $`a=b=x=(1C)/6`$, and $`r=(1+2C)/3`$, which can readily be seen to be a Werner state as in Eq. (24) or, equivalently,

$`\rho _{\mathrm{ii}}=\left(\begin{array}{cccc}\frac{2+C}{6}& 0& 0& \frac{1+2C}{6}\\ 0& \frac{1C}{6}& 0& 0\\ 0& 0& \frac{1C}{6}& 0\\ \frac{1+2C}{6}& 0& 0& \frac{2+C}{6}\end{array}\right).`$ (40)
Being the concurrence, $`C`$ is restricted to the interval $`[0,1]`$. The second solution is transcendental, but can be solved numerically.

In Fig. 8 we compare the four possible candidate solutions, and find that the global maximum is composed of only (i) and (ii).

We summarize the states at the frontier as follows:

$`\rho _{\mathrm{MEMS}:E_\mathrm{F},S_\mathrm{V}}=\{\begin{array}{cc}\rho _{\mathrm{ii}},\hfill & \mathrm{for}0CC^{};\hfill \\ \rho _\mathrm{i},\hfill & \mathrm{for}C^{}C1.\hfill \end{array}`$ (41)
Note the crossing point at $`(C,S_\mathrm{V})=(C^{},S_\mathrm{V}(C^{}))`$, at which extremality is exchanged, so the true frontier consists of two branches. It is readily seen that $`C^{}`$ is the solution of the equation $`S_\mathrm{V}(\rho _\mathrm{i}(C))=S_\mathrm{V}(\rho _{\mathrm{ii}}(C))`$, and the approximate numerical values of $`C^{}`$ and the corresponding $`S_\mathrm{V}^{}`$ are $`0.305`$ and $`0.741`$, respectively.

The resulting form of MEMS states is peculiar, in that, even at the crossing point of two branches on the entanglement-mixedness plane, the forms of matrices on the two branches are not equivalent (one is Rank 3, the other Rank 4). This is in contrast to the $`\{\mathrm{MEMS}:E_\mathrm{F},S_\mathrm{L}\}`$. This peculiarity can be partially understood from the plot of the two mixedness measures, Fig. 2: as the value of the von Neumann entropy rises, there are fewer and fewer Rank-3 entangled states, and above some threshold, no more Rank-3 states exist, let alone entangled Rank-3 states. There are, however, still entangled states of Rank 4. Hence, if Rank-3 states attain higher entanglement than Rank-4 states do when the entropy is low, a transition must occur between MEMS states of Rank 3 and Rank 4.

From Fig. 8 it is evident that beyond a certain value of the von Neumann entropy no entangled states exist. This value can be readily obtained by considering the MEMS state (40) at $`C=0`$,

$`\left(\begin{array}{cccc}\frac{1}{3}& 0& 0& \frac{1}{6}\\ 0& \frac{1}{6}& 0& 0\\ 0& 0& \frac{1}{6}& 0\\ \frac{1}{6}& 0& 0& \frac{1}{3}\end{array}\right),`$ (42)
for which $`S_\mathrm{V}=(1/2)\mathrm{log}_4(1/12)0.896`$.

As an aside, we mention a tantalizing but not yet fully developed analogy with thermodynamics Horodecki398 . In this analogy, one associates entanglement with energy and von Neumann entropy with entropy, and it is therefore tempting to regard the MEMS just derived as the analog of thermodynamic equilibrium states. If we apply the Jaynes principle to an ensemble in equilibrium with a given amount of entanglement then the most probable states are those MEMS shown above.

#### IV.3.2 Relative entropy of entanglement

Let us now find the frontier states for the case of relative entropy of entanglement. To do this, we first consider the Rank-3 states in Eq. (26), for which the relative entropy is given by Eq. (IV.2.2). For these states, the von Neumann entropy is given by

$$S_\mathrm{V}=\frac{1a+r}{2}\mathrm{log}\frac{1a+r}{2}a\mathrm{log}a\frac{1ar}{2}\mathrm{log}\frac{1ar}{2},$$
(43)
where the $`\mathrm{log}`$ function is taken to be base four. Even though the $`\mathrm{log}`$ functions in $`S_\mathrm{V}`$ and $`E_\mathrm{R}`$ use different bases, the stationary condition for the parameters $`r`$ and $`a`$ does not change, because the difference can be absorbed by a rescaling of the constraint-enforcing Lagrange multiplier. Thus, in maximizing $`E_\mathrm{R}`$ at fixed $`S_\mathrm{V}`$, we arrive at the stationarity condition

$$\mathrm{log}\frac{(1+a)^2r^2}{(1+a)^2}\mathrm{log}\frac{1ar}{1a+r}=\mathrm{log}\frac{1+a+r}{1+ar}\mathrm{log}\frac{(1a)^2r^2}{4a^2}.$$
(44)
We can solve for the parameter $`a`$ as a function of $`r[0,1]`$, at least numerically; the result is shown in Fig. 9, along with $`S_\mathrm{V}`$ and $`E_\mathrm{R}`$.

Turning to the Rank-4 case, it is straightforward, if tedious, to show that the Werner states satisfy the corresponding stationarity conditions. In order to ascertain which rank gives the MEMS for a given $`S_\mathrm{V}`$, we compare the stationary states of Rank 2, 3, and 4 in Fig. 10. Thus, we see that for $`S_\mathrm{V}S_\mathrm{V}^{}0.672`$ the frontier states are given by the Rank-3 states, whereas for $`S_\mathrm{V}S_\mathrm{V}^{}`$ the frontier states are given by the Werner states (24) with the parameter $`r`$ ranging from approximately 0.6059 down to 0. At the crossing point, $`(S_\mathrm{V}^{},E_\mathrm{R}^{})(0.672,0.124)`$, the MEMS undergo a discontinuous transition; recall that we encountered a similar phenomenon in the case of entanglement of formation vs. von Neumann entropy.

#### IV.3.3 Negativity

We saw in Sec. IV.2.3 that there is a pair of families of MEMS which differ in rank but give the identical frontier in the $`N`$ vs. $`S_\mathrm{L}`$ plane. It is interesting to see what happens for the combination of negativity and von Neumann entropy.

Once again, we begin with states of form (21), for which the negativity and the von Neumann entropy are given in Eqs. (32) and (37b), respectively. By making $`N`$ stationary at fixed $`S_\mathrm{V}`$, we are able to find only one solution (in addition to the Rank-2 candidate): $`a=b=x`$. Expressing the resulting density matrix, as we may, in terms of the single parameter $`r`$, we arrive at the following candidate for the frontier states:

$`\rho _{\mathrm{MEMS}:N,S_\mathrm{V}}=\left(\begin{array}{cccc}\frac{1+r}{4}& 0& 0& \frac{r}{2}\\ 0& \frac{1r}{4}& 0& 0\\ 0& 0& \frac{1r}{4}& 0\\ \frac{r}{2}& 0& 0& \frac{1+r}{4}\end{array}\right),`$ (45)
where $`0r1`$, i.e., the Werner states.

The resulting frontier in the negativity vs. von-Neumann-entropy plane is shown in Fig. 11 which, for comparison, also shows the curve for the Rank-2 candidate.

## V Concluding Remarks

In this Paper we have determined families of maximally entangled mixed states (MEMS, i.e., frontier states, which possess the maximum amount of entanglement for a given degree of mixedness). These states may be useful in quantum information processing in the presence of noise, as they have the maximum amount of entanglement possible for a given mixedness. We considered various measures of entanglement (entanglement of formation, relative entropy, and negativity) and mixedness (linear entropy and von Neumann entropy).

We found that the form of the MEMS depends heavily on the measures used. Certain classes of frontier states (such as those arising with either entanglement of formation or relative entropy of entanglement vs. the von Neumann entropy) behave discontinously at a specific point on the entanglement-mixedness frontier. Under most of the settings considered, we have been able to explicitly derive analytical forms for the frontier states.

For entanglement of formation and relative entropyand for most values of mixedness, we have found that the Rank-2 and Rank-3 MEMS have more entanglement than Werner states do. On the other hand, at fixed entropy no states have higher negativity than Werner states do. At small amounts of mixedness, the $`\{\mathrm{MEMS}:E_\mathrm{F},S_\mathrm{L}\}`$ states “lose” entanglement with increasing mixedness at a substantially lower rate than do the Werner states. However, when the entanglement is measured by the relative entropy, the difference in loss-rate is significantly smaller.

From Eq. (10) it is tempting to assert that in the case of LS entanglement vs. mixedness the frontier states should be the same as those in the case of entanglement of formation vs. mixedness. However, as we do not know whether Eq. (22) (up to local unitary transformations) exhausts all maximal states for LS entanglement, further investigation of this point is needed.

Having characterized the MEMS for various measures, it is worthwhile considering them from the perspective of Bell-inequality violations. To quantify the violation of Bell’s inequality, it is useful to consider the quantity

$$B\underset{\stackrel{}{a},\stackrel{}{a}^{},\stackrel{}{b},\stackrel{}{b}^{}}{\mathrm{max}}\left\{E(\stackrel{}{a},\stackrel{}{b})+E(\stackrel{}{a},\stackrel{}{b}^{})+E(\stackrel{}{a}^{},\stackrel{}{b})E(\stackrel{}{a}^{},\stackrel{}{b}^{})\right\},$$
(46)
where $`E(\stackrel{}{a},\stackrel{}{b})\stackrel{}{\sigma }\stackrel{}{a}\stackrel{}{\sigma }\stackrel{}{b}`$ and the vectors $`\stackrel{}{a}`$ and $`\stackrel{}{a}^{}`$ ($`\stackrel{}{b}`$ and $`\stackrel{}{b}^{}`$) are two different measuring apparatus settings for observer A (observer B). If $`B>2`$ then the corresponding state violates Bell’s inequality. For the density matrix of the form (21) it is straightforward Aravind95  to show that the quantity $`B`$ is given by

$`B=2\sqrt{\left(4(x+{\displaystyle \frac{r}{2}})1\right)^2+r^2}.`$ (47)
In Fig. 12, we plot $`B`$ vs. linear and von Neumann entropies for several families of frontier states. As a comparison, we also draw the violation by the following Rank-2 state (which is diagonal in the Bell basis):

$`\rho _\mathrm{B}r|\varphi ^+\varphi ^+|+(1r)|\varphi ^{}\varphi ^{}|,r[0,1].`$ (48)
This state, although not belonging to any of the families of frontier states derived previously, turns out to achieve the maximum possible violation, as a function of linear entropy. On the other hand, the Werner states appear to achieve maximal violation in the case of von Neumann entropy. Ekert’s application of Bell’s inequalities to quantum cryptography Ekert91 , together with the discussions of the present paragraph, suggests that MEMS may be relevant to quantum communication.

Another natural application for which entanglement is known to be a critical resource is quantum teleportation. How do these frontier MEMS teleport, compared with the Werner and Rank-2 Bell diagonal states? If we restrict our attention to high purity situations (i.e., to states with only a small amount of mixedness) then it is straightforward to show that, e.g., $`\{\mathrm{MEMS}:E_\mathrm{F},S_\mathrm{L}\}`$ states teleport average states better than the Werner states do, but worse than the Rank-2 Bell diagonal state does. Part of the explanation for this behavior is that standard teleportation is optimized for using Bell states as its core resource.

It is also interesting to note that for certain combinations of entanglement and mixedness measures, as well as the Bell inequality violation, the Rank-2 candidates fail to furnish MEMS. Thus, these states seem to be less useful than other MEMS. However, from the perspective of distillation, these states are exactly quasi-distillable Horodecki399 ; VerstraeteDehaeneDeMoor01 , and can be useful in the presence of noise because they can be easily distilled into Bell states.

## Acknowledgements

This material is based on work supported by NSF EIA01-21568 (TCW, PMG and PGK) and by the U.S. Department of Energy, Division of Material Sciences, under Award No. DEFG02-91ER45439, through the Frederick Seitz Materials Research Laboratory at the University of Illinois at Urbana-Champaign (TCW and PMG). KN and WJM acknowledge financial support from the European projects QUIPROCONE and EQUIP. PMG gratefully acknowledges the hospitality of the Unversity of Colorado at Boulder, where part of this work was carried out. TCW gratefully acknowledges the receipt of a Mavis Memorial Fund Scholarship.

## Appendix A Number of negative eigenvalues of the partial transpose of $`\rho `$

In this appendix we address the result that for $`C^2C^2`$ systems the partial transpose of any density matrix $`\rho `$ has at most one negative eigenvalue. In fact, we shall consider the result from two perspectives.

First, we build upon results (Theorem 3, in particular) contained in Ref. VerstraeteDehaeneDeMoor01 , from which it follows that it is sufficient to consider (i) Bell diagonal states and (ii) states of the form

$$\left(\begin{array}{cccc}a+c& 0& 0& d\\ 0& 0& 0& 0\\ 0& 0& bc& 0\\ d& 0& 0& ab\end{array}\right).$$
(49)
For the latter case, straightforward calculation shows that the partially transposed matrix can have one negative eigenvalue when $`d0`$ and that it does not have negative eigenvalues when $`d=0`$. For the former case, suppose that the Bell diagonal state has the four eigenvalues $`\lambda _1`$, $`\lambda _2`$, $`\lambda _3`$, $`\lambda _4`$, in nonincreasing order. Then it is straightforward to see that the corresponding partially transposed matrix has the four eigenvalues $`(\lambda _1+\lambda _2+\lambda _3\lambda _4)/2`$, $`(\lambda _1+\lambda _2\lambda _3+\lambda _4)/2`$, $`(\lambda _1\lambda _2+\lambda _3+\lambda _4)/2`$, $`(\lambda _1+\lambda _2+\lambda _3+\lambda _4)/2`$. Thus, it can have at most one negative eigenvalue.

A second perspective is provided by first invoking the LS decomposition (8) and then making a Schmidt decomposition, via the local unitary transformation $`U_AU_B`$, of the pure entangled part $`|\psi _e`$. In this way, the pure part becomes

$`U_AU_B|\psi _e=\sqrt{p}|00+\sqrt{1p}|11,`$ (50)
where $`1/2p1`$. Meanwhile, $`\rho `$ is transformed into

$`\rho ^{}`$ $`=`$ $`(U_AU_B)\rho (U_A^{}U_B^{})`$ (51)
$`=`$ $`\lambda \rho _s^{}+(1\lambda )\left(\begin{array}{cccc}p& 0& 0& \sqrt{p(1p)}\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ \sqrt{p(1p)}& 0& 0& 1p\end{array}\right),`$
where $`\rho _s^{}U_AU_B\rho _sU_A^{}U_B^{}`$ is still separable, and its partial transpose remains positive semi-definite.

Taking, then, the partial transpose, the density matrix is transformed into

$`\rho _{}^{}{}_{}{}^{\mathrm{T}_B}=\lambda \rho _{s}^{}{}_{}{}^{\mathrm{T}_B}+(1\lambda )\left(\begin{array}{cccc}p& 0& 0& 0\\ 0& 0& \sqrt{p(1p)}& 0\\ 0& \sqrt{p(1p)}& 0& 0\\ 0& 0& 0& 1p\end{array}\right),`$ (52)
and we note that the last matrix has eigenvalues $`\lambda _1=p`$, $`\lambda _2=\sqrt{p(1p)}`$, $`\lambda _3=1p`$, $`\lambda _4=\sqrt{p(1p)}`$ in (as $`p1/2`$) nonincreasing order. As the Hermitian matrix $`\rho _{s}^{}{}_{}{}^{\mathrm{T}_B}`$ retains positive-semi-definiteness, we can employ a well-known result in matrix analysis Horn\_Johnson  that for $`n\times n`$ Hermitian matrices $`A`$ and $`B`$, with $`B`$ being positive semi-definite, the eigenvalues of $`(A+B)`$ and $`A`$, when arranged in non-ascending order, obey

$`\lambda _k(A+B)\lambda _k(A)\text{for}k=1,2,\mathrm{},n.`$ (53)
Hence, identifying $`B`$ with $`\lambda \rho _{s}^{}{}_{}{}^{\mathrm{T}_B}`$ and $`A`$ with the product of $`(1\lambda )`$ and the matrix in Eq. (52), we immediately see that $`\lambda _1(\rho _{}^{}{}_{}{}^{\mathrm{T}_B})\lambda _2(\rho _{}^{}{}_{}{}^{\mathrm{T}_B})\lambda _3(\rho _{}^{}{}_{}{}^{\mathrm{T}_B})(1\lambda )(1p)0`$ and, thus, $`\rho _{}^{}{}_{}{}^{\mathrm{T}_B}`$ (or equivalently $`\rho ^{\mathrm{T}_B}`$) can have at most one negative eigenvalue. Thus, the negativity for $`C^2C^2`$ systems can then be written as $`N=2\mathrm{max}\{0,\lambda _4(\rho ^{\mathrm{T}_B})\}`$.


text:
# Local 4/5-Law and Energy Dissipation Anomaly in Turbulence

## 1 Introduction

Recently, Duchon and Robert  have established an energy balance relation for distributional solutions of the three-dimensional (3D) incompressible Euler equations. Their balance relation contains a “defect” or “anomaly” term, with an interesting connection to turbulence theory. Since the work of Duchon-Robert provides the point of departure of the present paper, it is appropriate to describe their theorems briefly here. For our purposes, there are three main results:

First, if $`𝐮L^3([0,T]\times 𝕋^3)`$ is a weak solution of the incompressible Euler equations on the 3-torus $`𝕋^3`$, then it is proved in  (Proposition 2) that the following local balance holds in the sense of distributions:

$$_t(\frac{1}{2}|𝐮|^2)+\mathbf{}\mathbf{}[(\frac{1}{2}|𝐮|^2+p)𝐮]=D(𝐮).$$
(1.1)
Here $`D(𝐮)`$ is a defect distribution which for classical solutions vanishes identically, implying local energy conservation. Duchon and Robert also establish various expressions for the defect term. In particular, they have shown that

$$D(𝐮)=𝒟\underset{\epsilon 0}{lim}\frac{1}{4}_{𝕋^3}\mathbf{}\phi _\epsilon (\mathbf{})\mathbf{}\delta 𝐮(\mathbf{})|\delta 𝐮(\mathbf{})|^2d^3\mathbf{}$$
(1.2)
where $`𝒟lim`$ means limit in the sense of distributions on $`[0,T]\times 𝕋^3`$, with $`\phi C_0^{\mathrm{}}(𝕋^3),`$ even, nonnegative with unit integral, $`\phi ^\epsilon (𝐱)=\epsilon ^3\phi (𝐱/\epsilon )`$, and $`\delta 𝐮(𝐱,t;\mathbf{})=𝐮(𝐱+\mathbf{},t)𝐮(𝐱,t)`$. This expression is remarkable because it is closely connected with an exact result in turbulence theory, the so-called “Kármán-Howarth-Monin relation”:

$$\mathbf{}\mathbf{}\delta 𝐮(\mathbf{})|\delta 𝐮(\mathbf{})|^2=4\overline{\epsilon }.$$
(1.3)
See , section 6.2.1. In this relation, $`\overline{\epsilon }=\nu |\mathbf{}𝐮^\nu |^2`$ is the mean energy dissipation for a Navier-Stokes solution $`𝐮^\nu `$, which is assumed to remain finite as viscosity $`\nu `$ tends to zero. We see that this result is essentially equivalent to the statement that the Duchon-Robert defect satisfy $`D=\overline{\epsilon }`$.

For a general distributional solution of Euler, there need be no connection of the defect with viscous dissipation, nor need it even be true that $`D(𝐮)0`$. However, a second theorem (Proposition 4) in  states that, if $`𝐮^\nu `$ is a Leray solution of incompressible Navier-Stokes equation for viscosity $`\nu `$ and if $`𝐮^\nu 𝐮`$ strong in $`L^3([0,T]\times 𝕋^3)`$ as $`\nu 0`$, then

$$D(𝐮)=𝒟\underset{\nu 0}{lim}[\nu |\mathbf{}𝐮^\nu |^2+D(𝐮^\nu )]$$
(1.4)
Since it is well-known that $`D(𝐮^\nu )0`$ for Leray solutions, thus also $`D(𝐮)`$ is a nonnegative distribution, i.e. a Radon measure (, Example 12.5). This theorem may be paraphrased as saying that strong viscosity solutions of the incompressible Euler equations are also dissipative solutions.

Finally, Duchon and Robert, under an additional hypothesis, establish an even simpler form of the defect distribution. With $`\omega `$ unit Haar measure on $`S^2`$, they define the function

$$S(𝐮,\mathrm{})(𝐱,t):=\frac{1}{\mathrm{}}_{S^2}𝑑\omega (\widehat{\mathbf{}})\delta u_L(𝐱,t;\mathbf{})|\delta 𝐮(𝐱,t;\mathbf{})|^2.$$
(1.5)
in $`L^1([0,T]\times 𝕋^3)`$. Here $`\delta u_L(𝐱,t;\mathbf{})=\widehat{\mathbf{}}\mathbf{}\delta 𝐮(𝐱,t;\mathbf{})`$ is the longitudinal velocity increment. Assuming that the following limit exists

$$S(𝐮)(𝐱,t):=𝒟\underset{\mathrm{}0}{lim}S(𝐮,\mathrm{})(𝐱,t)$$
(1.6)
Duchon and Robert , Section 5 show that

$$S(𝐮)=\frac{4}{3}D(𝐮).$$
(1.7)
This is a rigorous form of another well-known relation in turbulence theory, sometimes called the “4/3-law”:

$$\delta u_L(\mathbf{})|\delta 𝐮(\mathbf{})|^2\frac{4}{3}\overline{\epsilon }\mathrm{}.$$
(1.8)
See . It is well-known that the Kármán-Howarth-Monin relation reduces to the 4/3-law under conditions of local isotropy. This is achieved here by the angle average over the sphere in (1.5).

Such relations as (1.3) and (1.8) in turbulence theory go back to the original work of A. N. Kolmogorov . However, Kolmogorov in fact proved a relation involving only longitudinal velocity increments, the so-called “4/5-law”:

$$[\delta u_L(\mathbf{})]^3\frac{4}{5}\overline{\epsilon }\mathrm{}$$
(1.9)
This was established from the Navier-Stokes equations, under conditions of statistical homogeneity and local isotropy and with the assumption that energy dissipation remains finite in the zero-viscosity limit. It is our purpose here to establish an expression for the Duchon-Robert energy dissipation anomaly $`D(𝐮)`$ with exactly the form of Kolmogorov’s law. Our proof yields another expression which is related to the so-called “4/15-law”:

$$\delta u_L(\mathbf{})|\delta u_T(\mathbf{})|^2\frac{4}{15}\overline{\epsilon }\mathrm{}$$
(1.10)
Here $`\delta u_T(\mathbf{})=\widehat{𝒕}\mathbf{}\delta 𝐮(\mathbf{})`$ is a transverse velocity increment, with $`\widehat{𝒕}`$ any unit vector orthogonal to $`\widehat{\mathbf{}}`$. This relation is known to hold under the same conditions as the 4/5-law. The rigorous derivation of the 4/5- and 4/15-laws given here, under precisely stated assumptions, yields a result with a wider domain of validity than that of some previous rigorous derivations, such as that of Nie and Tanveer . In particular, the form of the 4/5-law established here—like the Duchon-Robert version of the 4/3-law—is local, in the sense that it relates third-order moments of velocity increments and viscous dissipation averaged over the same bounded spacetime region, not necessarily large.

In the following section we prove our main theorem. In a final discussion section we discuss briefly its physical significance and a compare it with previous results. Let us now make just a few remarks on notations: The symbol $`L^p`$ will be used below for $`L^p([0,T]\times 𝕋^3)`$. If $`F(𝐱,t)`$ is any spacetime distribution, we denote $`F^\epsilon (𝐱,t)=(\phi ^\epsilon F)(𝐱,t)`$, where $`\phi C_0^{\mathrm{}}(𝕋^3)`$ and $``$ is space convolution. Thus, $`F^\epsilon `$ remains, for each fixed $`𝐱`$, a distribution in $`t`$. Often below, as above, we omit the variables $`(𝐱,t)`$ where their presence is clear from the context.

## 2 The Main Theorem

We prove the following:

###### Theorem 1

Let $`𝐮L^3([0,T]\times 𝕋^3)`$ be a weak solution of the incompressible Euler equations on the 3-torus $`𝕋^3`$. Let $`\phi (\mathrm{})`$ be any $`C^{\mathrm{}}`$ function with compact support, nonnegative with unit integral, spherically symmetric, and let $`\phi ^\epsilon (\mathrm{})=\epsilon ^3\phi (\mathrm{}/\epsilon )`$. Finally, define longitudinal and transverse velocity increments as

$$\delta 𝐮_L(𝐱,t;\mathbf{})=(\widehat{\mathbf{}}\widehat{\mathbf{}})\delta 𝐮(𝐱,t;\mathbf{}),\delta 𝐮_T(𝐱,t;\mathbf{})=(\mathrm{𝟏}\widehat{\mathbf{}}\widehat{\mathbf{}})\delta 𝐮(𝐱,t;\mathbf{}),$$
(2.1)
Then, the following functions in $`L^1([0,T]\times 𝕋^3)`$

$$D_L^\epsilon (𝐮)=\frac{3}{4}_{𝕋^3}d^3\mathbf{}\left\{\mathbf{}\phi ^\epsilon (\mathrm{})\mathbf{}\delta 𝐮(\mathbf{})|\delta 𝐮_L(\mathbf{})|^2+\frac{2}{\mathrm{}}\phi ^\epsilon (\mathrm{})\delta u_L(\mathbf{})|\delta 𝐮_T(\mathbf{})|^2\right\},$$
(2.2)
and

$$D_T^\epsilon (𝐮)=\frac{3}{8}_{𝕋^3}d^3\mathbf{}\left\{\mathbf{}\phi ^\epsilon (\mathrm{})\mathbf{}\delta 𝐮(\mathbf{})|\delta 𝐮_T(\mathbf{})|^2\frac{2}{\mathrm{}}\phi ^\epsilon (\mathrm{})\delta u_L(\mathbf{})|\delta 𝐮_T(\mathbf{})|^2\right\},$$
(2.3)
both converge in the sense of distributions as $`\epsilon 0`$ to $`D(𝐮)`$, where the latter is the defect distribution in the local energy balance for $`𝐮`$,

$$_t(\frac{1}{2}|𝐮|^2)+\mathbf{}\mathbf{}[(\frac{1}{2}|𝐮|^2+p)𝐮]=D(𝐮),$$
(2.4)
established earlier by Duchon-Robert.

We shall prove this theorem in several steps.

The idea of the proof is to consider separate balance equations for the longitudinal and transverse components of the energy. We define first longitudinal and transverse velocities relative to a vector $`\mathbf{}`$:

$$𝐮_L(𝐱,t;\mathbf{})=(\widehat{\mathbf{}}\widehat{\mathbf{}})\mathbf{}𝐮(𝐱+\mathbf{},t),𝐮_T(𝐱,t;\mathbf{})=(\mathrm{𝟏}\widehat{\mathbf{}}\widehat{\mathbf{}})\mathbf{}𝐮(𝐱+\mathbf{},t).$$
(2.5)
Of course, $`𝐮_L+𝐮_T=𝐮`$. We define also mollified versions

$$𝐮_X^\epsilon (𝐱,t)=_{𝕋^3}d^3\mathbf{}\phi ^\epsilon (\mathrm{})𝐮_X(𝐱,t;\mathbf{}),X=L,T$$
(2.6)
It is easy to see that these satisfy the equations

$$_t𝐮_X^\epsilon +\mathbf{}\mathbf{}(𝐮𝐮_X)^\epsilon =\mathbf{}\mathbf{}𝚷_X^\epsilon ,X=L,T$$
(2.7)
distributionally in time, where

$$𝚷_L^\epsilon (𝐱,t)=_{𝕋^3}d^3\mathbf{}\phi ^\epsilon (\mathrm{})(\widehat{\mathbf{}}\widehat{\mathbf{}})p(𝐱+\mathbf{},t),𝚷_T^\epsilon (𝐱,t)=_{𝕋^3}d^3\mathbf{}\phi ^\epsilon (\mathrm{})(\mathrm{𝟏}\widehat{\mathbf{}}\widehat{\mathbf{}})p(𝐱+\mathbf{},t).$$
(2.8)
These equations can be simplified by the observation that

$$\mathbf{}\mathbf{}𝚷_X^\epsilon (𝐱,t)=\mathbf{}p_X^\epsilon (𝐱,t),X=L,T$$
(2.9)
where $`p_L^\epsilon ,p_T^\epsilon `$ are scalar functions defined by

$$p_X^\epsilon (𝐱,t)=_{𝕋^3}d^3\mathbf{}\phi _X^\epsilon (\mathrm{})p(𝐱+\mathbf{},t),X=L,T$$
(2.10)
with

$$\phi _L(\mathrm{})=\phi (\mathrm{})\phi _T(\mathrm{}),\phi _T(\mathrm{})=2_{\mathrm{}}^{\mathrm{}}\frac{\phi (\mathrm{}^{})}{\mathrm{}^{}}𝑑\mathrm{}^{}.$$
(2.11)
Note that $`\phi _L,\phi _T`$ are compactly supported and $`C^{\mathrm{}}`$ everywhere except at 0, where they have a mild (logarithmic) singularity. To prove (2.9), we use the elementary relation $`\mathbf{}\widehat{\mathbf{}}=\mathrm{𝟏}\widehat{\mathbf{}}\widehat{\mathbf{}}`$. A simple computation then gives, for example,

$$\mathbf{}\mathbf{}𝚷_L^\epsilon (𝐱,t)=_{𝕋^3}d^3\mathbf{}\left\{\frac{d\phi ^\epsilon }{d\mathrm{}}(\mathrm{})+\frac{2}{\mathrm{}}\phi ^\epsilon (\mathrm{})\right\}\widehat{\mathbf{}}p(𝐱+\mathbf{},t),$$
(2.12)
From its definition,

$$\mathbf{}\phi _L(\mathrm{})=\left\{\frac{d\phi ^\epsilon }{d\mathrm{}}(\mathrm{})+\frac{2}{\mathrm{}}\phi ^\epsilon (\mathrm{})\right\}\widehat{\mathbf{}}.$$
(2.13)
This gives easily $`\mathbf{}\mathbf{}𝚷_L^\epsilon =\mathbf{}p_L^\epsilon `$ in the sense of distributions. Because $`𝚷_L^\epsilon +𝚷_T^\epsilon =p^\epsilon \mathrm{𝟏}`$ and $`p_L^\epsilon +p_T^\epsilon =p^\epsilon `$, this yields also the relation $`\mathbf{}\mathbf{}𝚷_T^\epsilon =\mathbf{}p_T^\epsilon `$. Finally, we obtain the simpler equations

$$_t𝐮_X^\epsilon +\mathbf{}\mathbf{}(𝐮𝐮_X)^\epsilon =\mathbf{}p_X^\epsilon ,X=L,T$$
(2.14)
for $`𝐮_L^\epsilon `$ and $`𝐮_T^\epsilon `$.

We next observe that both $`𝐮_L^\epsilon `$ and $`𝐮_T^\epsilon `$ are divergence-free. In fact, a computation like that above shows that

$`\mathbf{}\mathbf{}𝐮_L^\epsilon (𝐱,t)`$ $`=`$ $`{\displaystyle _{𝕋^3}}d^3\mathbf{}\left\{{\displaystyle \frac{d\phi ^\epsilon }{d\mathrm{}}}(\mathrm{})+{\displaystyle \frac{2}{\mathrm{}}}\phi ^\epsilon (\mathrm{})\right\}\widehat{\mathbf{}}\mathbf{}𝐮(𝐱+\mathbf{},t),`$ (2.15)
$`=`$ $`{\displaystyle _{𝕋^3}}d^3\mathbf{}\mathbf{}\phi _L^\epsilon (\mathrm{})\mathbf{}𝐮(𝐱+\mathbf{},t).`$ (2.16)
If we integrate against any smooth test function $`\psi (𝐱)`$, we then get

$$_{𝕋^3}d^3𝐱\psi (𝐱)\mathbf{}\mathbf{}𝐮_L^\epsilon (𝐱,t)=_{𝕋^3}d^3\mathbf{}\phi _L^\epsilon (\mathrm{})\mathbf{}\mathbf{}(\psi 𝐮)(\mathbf{},t).$$
(2.17)
The latter is zero, since $`𝐮`$ is divergence-free (in the distributional sense). Thus, $`\mathbf{}\mathbf{}𝐮_L^\epsilon =0`$. Since $`𝐮_L^\epsilon +𝐮_T^\epsilon =𝐮^\epsilon ,`$ then also $`\mathbf{}\mathbf{}𝐮_T^\epsilon =0`$.

From the equations (2.14) for $`X=L,T`$, the incompressibility conditions for $`𝐮_L^\epsilon `$ and $`𝐮_T^\epsilon `$, and the Euler equations for $`𝐮`$ (in distribution sense), we derive the following balance equations

$$2_t(𝐮\mathbf{}𝐮_L^\epsilon )+\mathbf{}\mathbf{}[2(𝐮\mathbf{}𝐮_L^\epsilon )𝐮+((𝐮_L\mathbf{}𝐮_L)𝐮)^\epsilon (𝐮_L\mathbf{}𝐮_L)^\epsilon 𝐮+2p𝐮_L^\epsilon +2p_L^\epsilon 𝐮]=\frac{4}{3}D_L^\epsilon (𝐮),$$
(2.18)
and

$$2_t(𝐮\mathbf{}𝐮_T^\epsilon )+\mathbf{}\mathbf{}[2(𝐮\mathbf{}𝐮_T^\epsilon )𝐮+((𝐮_T\mathbf{}𝐮_T)𝐮)^\epsilon (𝐮_T\mathbf{}𝐮_T)^\epsilon 𝐮+2p𝐮_T^\epsilon +2p_T^\epsilon 𝐮]=\frac{8}{3}D_T^\epsilon (𝐮).$$
(2.19)
The basic identities used to derive these equations are

$`{\displaystyle _{𝕋^3}}d^3\mathbf{}\left\{\mathbf{}\phi ^\epsilon (\mathrm{})\mathbf{}\delta 𝐮(\mathbf{})|\delta 𝐮_L(\mathbf{})|^2+{\displaystyle \frac{2}{\mathrm{}}}\phi ^\epsilon (\mathrm{})\delta u_L(\mathbf{})|\delta 𝐮_T(\mathbf{})|^2\right\}`$ (2.20)
$`={\displaystyle _{𝕋^3}}d^3\mathbf{}{\displaystyle \frac{}{\mathrm{}_k}}\{\widehat{\mathrm{}}_i\widehat{\mathrm{}}_j\phi ^\epsilon (\mathrm{})\}\delta u_i(\mathbf{})\delta u_j(\mathbf{})\delta u_k(\mathbf{})`$ (2.21)
$`={\displaystyle \frac{}{x_k}}[((𝐮_L\mathbf{}𝐮_L)u_k)^\epsilon (𝐮_L\mathbf{}𝐮_L)^\epsilon u_k]+2u_i{\displaystyle \frac{}{x_k}}[((u_{Li}u_k)^\epsilon u_{Li}^\epsilon u_k],`$ (2.22)
and

$`{\displaystyle _{𝕋^3}}d^3\mathbf{}\left\{\mathbf{}\phi ^\epsilon (\mathrm{})\mathbf{}\delta 𝐮(\mathbf{})|\delta 𝐮_T(\mathbf{})|^2{\displaystyle \frac{2}{\mathrm{}}}\phi ^\epsilon (\mathrm{})\delta u_L(\mathbf{})|\delta 𝐮_T(\mathbf{})|^2\right\}`$ (2.23)
$`={\displaystyle _{𝕋^3}}d^3\mathbf{}{\displaystyle \frac{}{\mathrm{}_k}}\{(\delta _{ij}\widehat{\mathrm{}}_i\widehat{\mathrm{}}_j)\phi ^\epsilon (\mathrm{})\}\delta u_i(\mathbf{})\delta u_j(\mathbf{})\delta u_k(\mathbf{})`$ (2.24)
$`={\displaystyle \frac{}{x_k}}[((𝐮_T\mathbf{}𝐮_T)u_k)^\epsilon (𝐮_T\mathbf{}𝐮_T)^\epsilon u_k]+2u_i{\displaystyle \frac{}{x_k}}[((u_{Ti}u_k)^\epsilon u_{Ti}^\epsilon u_k].`$ (2.25)
Now we study the convergence of the lefthand sides of (2.18) and (2.19) as $`\epsilon 0`$.

First we show that $`𝐮_L^\epsilon \frac{1}{3}𝐮`$ and $`𝐮_T^\epsilon \frac{2}{3}𝐮`$ strong in $`L^3`$ as $`\epsilon 0`$. Because of the spherical symmetry, compact support and unit integral of $`\phi ^\epsilon `$,

$$_{𝕋^3}d^3\mathbf{}\phi ^\epsilon (\mathrm{})\widehat{\mathbf{}}\widehat{\mathbf{}}=\frac{1}{3}\mathrm{𝟏},$$
(2.26)
for sufficiently small $`\epsilon `$. From this result and the definition of $`𝐮_L^\epsilon `$ it follows that

$$𝐮_L^\epsilon (𝐱,t)\frac{1}{3}𝐮(𝐱,t)=_{𝕋^3}d^3\mathbf{}\phi ^\epsilon (\mathrm{})\widehat{\mathbf{}}\widehat{\mathbf{}}\mathbf{}[𝐮(𝐱+\mathbf{},t)𝐮(𝐱,t)].$$
(2.27)
Thus,

$$𝐮_L^\epsilon \frac{1}{3}𝐮_{L^3}_{𝕋^3}d^3\mathbf{}\phi ^\epsilon (\mathrm{})𝐮(+\mathbf{})𝐮_{L^3}.$$
(2.28)
Because $`𝐮L^3`$, it follows by a standard approximation argument that $`𝐮(+\mathbf{})𝐮_{L^3}0`$ as $`\mathrm{}0`$. Hence, it follows from (2.28) that $`𝐮_L^\epsilon \frac{1}{3}𝐮_{L^3}0`$ as $`\epsilon 0`$, as was claimed. Because $`𝐮_T^\epsilon =𝐮^\epsilon 𝐮_L^\epsilon `$, we also deduce that $`𝐮_T^\epsilon \frac{2}{3}𝐮_{L^3}0`$.

Next we show that $`p_L^\epsilon \frac{1}{3}p`$ and $`p_T^\epsilon \frac{2}{3}p`$ strong in $`L^{3/2}`$ as $`\epsilon 0`$. First, observe that $`p=()^1_i_j(u_iu_j)`$, so that $`pL^{3/2}`$ by the Calderón-Zygmund inequality. From the definitions of $`\varphi _T`$ and $`p_T^\epsilon `$ it follows easily that

$$p_T^\epsilon (𝐱,t)=\frac{2}{3}4\pi _0^{\mathrm{}}\mathrm{}^2𝑑\mathrm{}^{}\phi ^\epsilon (\mathrm{}^{})\frac{1}{|B_{\mathrm{}^{}}|}_B_{\mathrm{}^{}}d^3\mathbf{}p(𝐱+\mathbf{},t),$$
(2.29)
where $`B_{\mathrm{}}`$ is the ball of radius $`\mathrm{}`$ centered at the origin and $`|B_{\mathrm{}}|`$ is its volume. Because $`\phi ^\epsilon `$ has unit integral,

$$p_T^\epsilon (𝐱,t)\frac{2}{3}p(𝐱,t)=\frac{2}{3}4\pi _0^{\mathrm{}}\mathrm{}^2𝑑\mathrm{}^{}\phi ^\epsilon (\mathrm{}^{})\frac{1}{|B_{\mathrm{}^{}}|}_B_{\mathrm{}^{}}d^3\mathbf{}[p(𝐱+\mathbf{},t)p(𝐱,t)],$$
(2.30)
and thus

$$p_T^\epsilon \frac{2}{3}p_{L^{3/2}}\frac{2}{3}4\pi _0^{\mathrm{}}\mathrm{}^2d\mathrm{}^{}\phi ^\epsilon (\mathrm{}^{})\frac{1}{|B_{\mathrm{}^{}}|}_B_{\mathrm{}^{}}d^3\mathbf{}p(+\mathbf{})p_{L^{3/2}}.$$
(2.31)
Consequently, $`p_T^\epsilon \frac{2}{3}p_{L^{3/2}}0`$ as $`\epsilon 0`$, as was to be proved. Since $`p_L^\epsilon =p^\epsilon p_T^\epsilon `$, also $`p_L^\epsilon \frac{1}{3}p_{L^{3/2}}0`$.

Entirely similar arguments show that $`((𝐮_X\mathbf{}𝐮_X)𝐮)^\epsilon (𝐮_X\mathbf{}𝐮_X)^\epsilon 𝐮0`$ strong in $`L^1`$ as $`\epsilon 0`$, for $`X=L,T.`$ We leave that to the reader.

Finally we see that the the lefthand side of equation (2.18) converges as $`\epsilon 0`$ to the quantity

$$\frac{4}{3}\left\{_t(\frac{1}{2}|𝐮|^2)+\mathbf{}\mathbf{}\left[(\frac{1}{2}|𝐮|^2+p)𝐮\right]\right\}$$

in the sense of distributions, and that the lefthand side of (2.19) converges to $`\frac{8}{3}`$ times the same quantity in the curly brackets, also in the distribution sense. However, by the result of Duchon-Robert , that quantity in the brackets is equal to minus the defect distribution $`D(𝐮)`$. We conclude then that

$$D_X^\epsilon (𝐮)D(𝐮)$$
(2.32)
as $`\epsilon 0`$ for both $`X=L,T`$, in the sense of distributions. $`\mathrm{}`$

We now prove the following:

###### Corollary 1

Assume that the functions $`S_L(𝐮,\mathrm{}),S_T(𝐮,\mathrm{})L^1([0,T]\times 𝕋^3)`$ defined by

$$S_L(𝐮,\mathrm{})=\frac{1}{\mathrm{}}_{S^2}𝑑\omega (\widehat{\mathbf{}})[\delta u_L(\mathbf{})]^3$$
(2.33)
and

$$S_T(𝐮,\mathrm{})=\frac{1}{\mathrm{}}_{S^2}𝑑\omega (\widehat{\mathbf{}})\delta u_L(\mathbf{})|\delta 𝐮_T(\mathrm{})|^2$$
(2.34)
have limits

$$S_X(𝐮)=𝒟\underset{\mathrm{}0}{lim}S_X(𝐮,\mathrm{}),X=L,T.$$
(2.35)
Then,

$$S_L(𝐮)=\frac{4}{5}D(𝐮),S_T(𝐮)=\frac{8}{15}D(𝐮).$$
(2.36)
The proof is quite straightforward. An easy computation for a spherically symmetric test function gives

$$\frac{4}{3}D_L^\epsilon (𝐮)=4\pi _0^{\mathrm{}}𝑑\mathrm{}\left[\mathrm{}^3\phi ^{}(\mathrm{})S_L(𝐮,\epsilon \mathrm{})+2\mathrm{}^2\phi (\mathrm{})S_T(𝐮,\epsilon \mathrm{})\right]$$
(2.37)
and

$$\frac{8}{3}D_T^\epsilon (𝐮)=4\pi _0^{\mathrm{}}𝑑\mathrm{}\left[\mathrm{}^3\phi ^{}(\mathrm{})2\mathrm{}^2\phi (\mathrm{})\right]S_T(𝐮,\epsilon \mathrm{}).$$
(2.38)
Taking the limit as $`\epsilon 0`$, using Theorem 1, the normalization of the test function, and the above hypotheses, we see that

$$\frac{4}{3}D(𝐮)=3S_L(𝐮)+2S_T(𝐮)$$
(2.39)
and

$$\frac{8}{3}D(𝐮)=5S_T(𝐮).$$
(2.40)
Solving this linear system gives the result. $`\mathrm{}`$

In the definition of $`S_T(𝐮,\mathrm{})`$ in (2.34) we could instead have replaced $`|\delta 𝐮_T|^2`$ by the square magnitude of a transverse component $`\delta u_T=\widehat{𝒕}\mathbf{}\delta 𝐮`$, where $`\widehat{𝒕}`$ is any unit vector perpendicular to $`\widehat{\mathbf{}}`$. Call such a quantity $`\stackrel{~}{S}_T(𝐮,\mathrm{})`$. However, $`|\delta 𝐮_T|^2=|\delta u_T|^2+|\delta u_T^{}|^2,`$ where $`\delta u_T^{}=(\widehat{𝒕}\mathbf{\times }\widehat{\mathbf{}})\mathbf{}\delta 𝐮.`$ Since both transverse components give an equal contribution in the spherical average, $`S_T(𝐮,\mathrm{})=2\stackrel{~}{S}_T(𝐮,\mathrm{}).`$ Under the assumptions of the corollary, $`\stackrel{~}{S}_T(𝐮,\mathrm{})`$ then also has a limit distribution $`\stackrel{~}{S}_T(𝐮)`$ as $`\mathrm{}0`$, and

$$\stackrel{~}{S}_T(𝐮)=\frac{4}{15}D(𝐮).$$
(2.41)
This is the conventional statement of the “4/15-law”.

## 3 Discussion

If we couple our Corollary 1 with Duchon-Robert’s characterization of $`D(𝐮)`$ by the inviscid limit, , Proposition 4, then we arrive at essentially the following statement: For any “nice” spacetime region $`R𝕋^3\times [0,T]`$, let

$$\epsilon _R=\underset{\nu 0}{lim}\frac{1}{|R|}_Rd^3𝐱𝑑t\epsilon (𝐱,t)$$
(3.1)
with $`\epsilon (𝐱,t):=\nu |\mathbf{}𝐮^\nu |^2`$ (assuming for simplicity that $`D(𝐮^\nu )=0`$). If first $`\nu 0`$ and next $`\mathrm{}0`$, then

$$[\delta u_L(\mathbf{})]^3_{ang,R}\frac{4}{5}\epsilon _R\mathrm{}$$
(3.2)
with $`_{ang,R}`$ denoting an average of $`(\widehat{\mathbf{}},𝐱,t)`$ over $`S^2\times R`$. To be precise, the relation should hold in the sense of averaging in $`𝐱,t`$ against smooth test functions $`\psi (𝐱,t)`$ rather than against a characteristic function $`1_R(𝐱,t)`$. The latter sense is stronger, because, if the result holds for characteristic functions of Borel sets $`R`$, then, by approximation arguments, it holds also for smooth functions. In fact, by the portmanteau theorem , it is enough to show that it holds for all bounded Borel sets $`R`$ such that $`D(R)=0`$, i.e. with no energy dissipation concentrating on the boundary. We see from these remarks that the meaning of the Corollary 1 is that the scaling observed in the Kolmogorov “4/5-law”  should hold with just local averaging in spacetime and an angular average over the direction of the separation vector. Note that, if the solution of the Navier-Stokes equations at finite viscosity is not regular, then the same relation would hold with the modified definition $`\epsilon (𝐱,t):=\nu |\mathbf{}𝐮^\nu |^2+D(𝐮^\nu )`$, including the dissipation from the singularities.

Another recent work by Nie and Tanveer  has given a rigorous derivation of both the 4/5- and 4/3-laws under conditions similar to those of Duchon-Robert and of the present paper. The proof of those laws by Nie and Tanveer is for individual solutions of the Navier-Stokes equations (assumed to be strong) without any statistical averaging, and with no assumptions of homogeneity or isotropy. In , as here, the averaging is over spacetime and over directions of the separation vector. However, Nie and Tanveer considered the space average over the entire domain and time average over an interval $`[0,T]`$ with $`T\mathrm{}`$. In contrast, the present result holds for any spacetime domain, of arbitrary extent, so long as the Reynolds number is sufficiently high, and also for solutions of the Navier-Stokes equations that may be singular. In this respect, our result is stronger than that of . On the other hand, our proof and that in  give no indication how large the Reynolds number must be taken to approach the limit, whereas Nie and Tanveer establish precise error bounds for their result.

There is a slight resemblance of our and Duchon-Robert’s local results with the “refined similarity hypothesis”that Kolmogorov postulated in 1962 . Kolmogorov considered averages $`\epsilon _B`$ of dissipation over balls $`B𝕋^3`$ and the structure functions obtained by averaging velocity increments $`\delta 𝐮(𝐱,t;\mathbf{})`$ conditioned upon the local dissipation averaged over a ball $`B_{\mathrm{}}(𝐱)`$ of diameter $`\mathrm{}`$ centered at the midpoint $`\frac{1}{2}(𝐱+\mathbf{})`$. He postulated that these might be calculated by his original 1941 theory, e.g. for longitudinal differences. Hence, he obtained for unconditioned ensemble averages

$$[\delta u_L(\mathbf{})]^{3n}C_n\epsilon _B_{\mathrm{}}^n\mathrm{}^n$$
(3.3)
Our scaling relation (3.2) is supposed to hold pathwise, for individual realizations. Hence, if it holds—along with a uniform and integrable bound for $`\mathrm{}0`$— then one can infer that

$$\left[(\delta u_L(\mathbf{}))^3_{ang,R}\right]^n\left(\frac{4}{3}\right)^n\epsilon _R^n\mathrm{}^n.$$
(3.4)
The similarity is apparent. Aside from the fact that we must integrate over time as well as space, the major and significant difference is that in the relation (3.4) the region $`R`$ must be fixed as $`\mathrm{}0`$, while in Kolmogorov’s hypothesis the ball $`B_{\mathrm{}}`$ itself shrinks to zero in the limit. Our results and that of Duchon-Robert say nothing about intermittency.

Although we have proved a rigorous theorem, it is only established under various hypotheses. We regard those as plausible, but it is nevertheless of some interest to inquire about the feasibility of numerical and/or experimental tests of the local results proved here and in . In , the 4/5-law has been verified both in a numerical simulation via a space-average over the domain and also in experimental hot-wire data. However, in the same simulation, the 4/15-law is not satisfied (S. Chen, private communication). Similar violations of the 4/15-law were observed in experimental X-wire data from an atmospheric boundary layer . In both cases, anisotropy of the flow may be suspected of vitiating the result. In that case, angle-averaging as employed in the rigorous theorems may improve the agreement. On the other hand, a more recent numerical study  with a $`1024^3`$ simulation has found that the regimes of validity of both the 4/5- and 4/3-laws, in even their global form, are established only very slowly as the Reynolds number is increased. Such a slow approach to asymptotia makes a direct test of the local results, especially a verification of the numerical prefactor, rather more difficult.

Acknowledgements: This work, LAUR 02-4761, was carried out in part at the Los Alamos National Laboratory under the auspices of the Department of Energy and supported by LDRD-ER 2000047. The work was begun at the Erwin Schroedinger International Institute for Mathematical Physics during its Program on Developed Turbulence in spring 2002. I wish to acknowledge several of my colleagues there for fruitful discussions, in particular Jean Duchon, Susan Kurien, and Christos Vassilicos. I would also like to thank Mark Taylor of Los Alamos National Laboratory and Shiyi Chen at Johns Hopkins University for information about their numerical simulations.


text:
# Does Betelgeuse have a Magnetic Field?

## 1. Introduction

While the cool super-giant star Betelgeuse ($`\alpha `$ Orionis) is among the stars with the largest apparent diameters—corresponding to a radius in the range 600–800 $`\mathrm{R}_{}`$—fundamental stellar parameters for this red M1–2 Ia–Iab star are by no means well-defined. Recently Freytag and collaborators (e.g. Freytag, Steffen, & Dorch 2002) performed detailed numerical three-dimensional radiation-hydrodynamic simulations of the outer convective envelope and atmosphere of the star under realistic physical assumptions. They try to determine if its observed brightness variations may be understood by convective motions within the star’s atmosphere: the resulting models are largely successful in explaining the observations as a consequence of giant-cell convection on the stellar surface, very dissimilar to solar convection. These detailed simulations bring forth the possibility of solving another question regarding the nature of Betelgeuse, and of super-giants in general; namely whether these stars may harbor magnetic activity, which in turn may also contribute to their variability. A possible astrophysical dynamo in Betelgeuse would most likely be very different from those thought to operate in solar type stars, both due to its slow rotation, and to the fact that only a few convection cells are present at its surface at any one time.

## 2. Dynamo model

The basic ansatz of the approach in this paper is that the input prescribed flow-field is taken at face value; i.e. that the velocity ceteris paribus represents the true quantity in the real Betelgeuse.

We solve the induction equation for the prescribed velocity field, i.e. within the kinematic MHD approximation, which is valid when the magnetic field is weak:

$$\frac{𝐁}{t}=\times (𝐮\times 𝐁)+\eta ^2𝐁,$$
(1)
where $`𝐁`$ is the magnetic field (flux density), $`𝐮`$ is the prescribed velocity field, and $`\eta `$ is the magnetic diffusivity (the resistivity).

When $`𝐁`$ becomes comparable to the equipartition value of the convective flows B$`{}_{\mathrm{eq}}{}^{}=\mathrm{u}\sqrt{\mu \rho }`$, non-linear effects becomes important through the back-reaction of the Lorentz force on the flow. We are not able to model this non-linear behavior exactly within the present approach—the flow field is taken from a long since done calculation—instead we attempt an ad hoc strategy: as a first step towards including non-linearity, in a few cases we replace the flow $`𝐮`$ in Eq. (1) by a flow $`𝐮_\mathrm{q}`$ quenched by the magnetic field:

$$𝐮_\mathrm{q}=𝐮\mathrm{exp}\alpha (e_\mathrm{m}/e_\mathrm{k})^2,$$
(2)
where $`\alpha `$ is a constant, $`e_\mathrm{m}`$ the magnetic energy density, and $`e_\mathrm{k}`$ is the average kinetic energy density, taken to be constant in both space and time.

This quenching thus reduces the velocity amplitude at the locations where the magnetic energy becomes comparable to the kinetic energy of the fluid flow. Hence it reduces the growth rate of the magnetic field in these regions, causing the total magnetic energy E$`{}_{\mathrm{m}}{}^{}=_\mathrm{V}\mathrm{e}_\mathrm{m}\mathrm{dV}`$ to saturate. We have chosen $`\alpha =1.75`$, which ensures the flow to be effectively unquenched until field strengths of $`0.1\mathrm{B}_{\mathrm{eq}}`$, while the flow amplitude completely vanishes at field strengths of $`1.2\mathrm{B}_{\mathrm{eq}}`$. Note that the constant $`e_\mathrm{k}`$ merely introduces a scaling factor in the solution to $`𝐁`$, and that its actual value is irrelevant to the present formulation of the problem. In future work the simple quenching Eq. (2) will be replaced with an expression taking into account the geometries of the field and flow.

### 2.1. Numerical method

The employed numerical scheme for the kinematic dynamo simulations is based on the staggered grid method by Galsgaard, Nordlund and others (e.g. Galsgaard & Nordlund 1997): it uses 6th order staggering operators, 5th order centering routines, and a 3rd order Hyman predictor-corrector time-stepping. The same code has been used in the context of dynamo action previously, to study both kinematic dynamo action (Dorch 2000), as well as non-linear turbulent dynamos (Archontis, Dorch, & Nordlund 2002), and recently to study kinematic dynamo action by convective flows in M-type dwarf stars (Dorch & Ludwig 2002).

We implement a perpendicular magnetic field boundary condition by introducing symmetry conditions on the magnetic field across ghost zones at the boundaries in all three dimensions. One could argue that ideally radial or potential field boundary condition are more physical, and it is indeed the plan to use such realistic boundary conditions in future work. However, we have tested different simple boundaries and the results in terms of integral quantities such as e.g. magnetic energy seem quite robust.

Likewise, we have employed different types of initial conditions for the magnetic field and find that after a short transient, the behavior of the system does not depend on the particular choice of topology, as it will become apparent from the subsequent discussion of our results. The simulations are initiated with either a unidirectional weak field, or a periodic weak field with a large number of null-points (an ABC-like topology, see e.g. Dorch 2000).

### 2.2. Flow field input data

The input flow field results from the star-in-a-box models of co-author Freytag and collaborators (see e.g. Freytag, Steffen, & Dorch 2002). The particular dataset used here is from a rather coarse model with $`127^3`$ grid points; the model identification designation is dst33gm06n03.

The giant convection cells are so large that only a few cells are present at the surface at the same time; hence the situation is very different from that of dwarf stars such as the Sun, where there can be thousands of cells present at the surface. Betelgeuse is only slowly rotating and the star-in-a-box model does not include rotational forces: it follows that the flows are not very helical. In fact the relative kinetic helicity $`\omega 𝐮/(\omega _{\mathrm{rms}}\mathrm{u}_{\mathrm{rms}})`$ is only on the order of 0.02 (where $`\omega =\times 𝐮`$ is the vorticity).

The flow field data comprises 120 snapshots covering 7.5 years of giant cell convection. This is actually not the full duration of that particular simulation, but we had to limit the amount of input data, due to lack of available disk space. The time step resulting from the Hyman predictor-corrector scheme used when solving Eq. (1) is typically 25 times smaller than the interval between the flow field snapshots; interpolation at each time step is hence necessary and sufficiently smooth behavior is achieved using a simple first order interpolation routine.

To be able to study longer time sequences than the 7.5 year extent of the input data, we cyclically re-use the input flow, effectively introducing a 7.5 year periodicity in the flow; while such a periodicity arguably is not observed, nor likely in Betelgeuse, it allows us to study long term effects related to diffusion on global scales.

### 2.3. Diffusion and magnetic Reynolds number

Dynamo action by flows are often studied in the limit of increasingly large magnetic Reynolds numbers Re$`{}_{\mathrm{m}}{}^{}=\mathrm{}\mathrm{U}/\eta `$, where $`\mathrm{}`$ and U are characteristic length and velocity scales. Most astrophysical systems are highly conducting (yielding small magnetic diffusivities/resistivities $`\eta `$), and their dimensions are huge; consequently values of Rem are huge too. Seemingly odd exceptions are e.g. cool M-type dwarf stars (see Dorch & Ludwig 2002) that have atmospheres a hundred times more neutral than the Sun. Betelgeuse is not an exception however; most parts of the star is better conducting than the solar surface layers, which has a magnetic diffusivity of the order of $`\eta 10^4`$ m2/s.

Figure 1 shows the average Spitzer’s resistivity as a function of radius in the model of Betelgeuse: Spitzer’s formula (e.g. Schrijver & Zwaan 2000) assumes complete ionization and hence the precise values of $`\eta `$ are uncertain in the outer parts of the star, where the atmosphere borders on neutral. There is some uncertainty connected also with defining the most important length scale of the system, but taking $`\mathrm{}`$ to be 10% of the radial distance R from the center (a typical scale of the giant cells), and U $`=\mathrm{u}_{\mathrm{RMS}}`$ along the radial direction yields Re$`{}_{\mathrm{m}}{}^{}=10^{10}`$$`10^{12}`$ in the interior part of the star where R $`700`$ R.

Our numerical approach invokes quenching of the magnetic diffusivity by the convergence of the flow-field across the magnetic field, so that the magnetic Reynolds number Rem is large in the bulk of the flow (this diffusion quenching should not be confused with the quenching of the flow given by Eq. 2). The primary advantage of this approach is that diffusion is confined to the small regions where it is in fact needed, in order to resolve the smallest magnetic structures on the numerical grid (where the magnetic field gradients are large): typically we set the minimum value of Rem larger than a few hundred.

## 3. Results

There seems to be some disagreement as to what one should require for a system to be an astrophysical dynamo. Several ingredients can be considered to be necessary in order for a system to be a “true” dynamo; we believe that the following four should suffice.

1. The flows must stretch, twist and fold the magnetic field lines.
2. Reconnection must take place to render the above processes irreversible (i.e. diffusion is needed locally).
3. The weak magnetic field must be circulated to the locations where flow can do (Lorentz) work upon it.
4. The total volume magnetic energy Em must increase (if a kinematic dynamo).

The conjecture that these necessary ingredients are sufficient is based largely on our experience from idealized kinematic and non-linear dynamo models (Archontis & Dorch 1999; Dorch 2000; Archontis, Dorch, & Nordlund 2002).

The present contribution presents however only a preliminary study of a possible Betelgeusian dynamo, and we shall be dealing mainly with the last of these four ingredients: namely the question of exponential growth of Em. Identification of the mode of operation of the dynamo will be presented elsewhere.

In general we obtain dynamo action when the specified minimum value of Rem is larger than approximately 500: at lower values of Rem the total magnetic energy decays throughout. In the following we report on results for cases with $`min(\mathrm{Re}_\mathrm{m})=800`$. In that case Rem is much larger (on the order of $`10^{10}`$) in the bulk of the flow.

Figure 2 shows the evolution of the total volume magnetic energy Em during the 7.5 year time interval covered by the flow field data: initially there is a short increase in Em during the first half year, which is followed by a two year decline, and a return to exponential growth in the reminder of the simulation. One could worry that the particular choice of magnetic field initial condition would influence the results, but as already hinted this is not the case, as is evident in Figure 3, which shows the evolution of Em over 90 years. In this model, the flow field data was re-used 12 times (patched together at the ends by interpolated). The first 4 recyclings happens to begin with different magnetic field configurations, and they can thus be thought of as corresponding to simulations with different initial conditions: in all cases there is an over-all exponential growth of the total magnetic energy. In the last 9 cycles until the simulation was terminated, the growth of the magnetic field proceeds in a very regular (artificial, one might add) manner, with a well defined average exponential growth rate corresponding to a growth time of 14.7 years.

In Figure 2 we have overplotted the average evolution of Em in the last 9 cycles of the simulation: in these cycles the behavior of Em is completely regular, and we speculate the it corresponds to some eigen-mode of the dynamo (determined by Eq. 2), i.e. the field geometry at the beginning of each cycle is the one that belongs to the largest growth rate, which at $`min(\mathrm{Re}_\mathrm{m})=800`$ turns out to be $`0.07`$ yr-1.

No exponential growth can go on forever and eventually the magnetic energy amplification must saturate. To model this non-linear effect, we ran a simulation identical to the kinematic model yielding the result in Figure 3, but including the flow quenching of Eq. (2). Figure 4 shows the evolution of a saturating dynamo, after an initial 25 year exponential growth (identical to that of the purely kinematic simulation): in that case the flow field cycles become less apparent in the evolution of Em, and a additional long term variations become visible—the dominant period seems to be about 35 years (more than 4.5 cycles). The simulation was terminated after 300 years in Betelgeusian time. One may worry that we have not run the simulation for long enough compared to some dominate diffusion time scale. If diffusion takes place on a time scale related to the largest scale of the system (see e.g. Brandenburg 2001), then the characteristic diffusion time would be $`\tau _\mathrm{d}=\mathrm{R}^2/4\pi ^2\eta 210^{10}`$ years, i.e. longer than the age of the Universe. However, scaling this with our minimum Re$`{}_{\mathrm{m}}{}^{}=800`$ (to its actual value that is on the order of $`10^{10}`$) yields a diffusion time of $`\tau _\mathrm{d}1500`$ years: that time scale does not presently seem completely inconsistent with our results (see Figure 4), but the simulation covers only a fraction of $`\tau _\mathrm{d}`$.

The level of saturation of Em is effectively set by the choice of the constant kinetic energy density in Eq. (2): thus the model does not contribute any additional knowledge about which field strengths to expect in Betelgeuse.

But we can, however, study the geometry of the magnetic field that the dynamo generates. Figure 5 (left) shows a volume rendering of isosurfaces of magnetic field structures with a high field strength relative to the maximum: the field becomes concentrated into elongated structures much thinner than the scale of the giant convection cells, but perhaps due to the very dynamic nature of the convective flows, no “intergranular network” is formed (to use solar terminology). The field is highly intermittent (see PDF in Figure 6, left), i.e. only a small fraction of the volume carries the strongest structures. The average magnetic field distribution is illustrated in Figure 5 (right), where it is evident that the field is stratified and decreases in strength from the center of the star to a rather sharp “magnetic surface” at R $`750`$ R. We speculate that the stratification results from the magnetic pumping effect working on the weakest part of the field (cf. Dorch & Nordlund 2001).

The magnetic structures are well resolved, with maximum power on the largest scales above 200 R (Figure 6, right). Additionally we observe a slight trend in the topology of the field; fields near (but still below) the surface of the star are predominantly horizontally aligned, while those in deeper layers are radial.

## 4. Concluding remarks

Based on the results presented here, we may not say conclusively if Betelgeuse does have a magnetic field, of course. The results are tentative and should be used with caution. But we may say that it seems that it can indeed have a presently unobserved magnetic field. The dynamo of Betelgeuse may be characterized as belonging to the class called “local small-scale dynamos” another example of which is the proposed dynamo action in the solar photosphere that may possibly be responsible for the formation of small-scale flux tubes (cf. Cattaneo 1999, but also the discussion by Stein, ibid). However, in the case of Betelgeuse this designation is less meaningful since the generated magnetic field is both global and large-scale.

The future developments of this project will involve firstly using longer time sequences of the input flow field, to avoid having to rely on recycling. Secondly, simulations with higher numerical resolution is currently being carried out (see Freytag & Finnsson, 2002, and Freytag & Mizuno-Wiedner, 2002), which will allow larger runs with higher magnetic Reynolds numbers (i.e. smaller magnetic structures can form). Thirdly, it will be a priority to apply a more realistic quenching expression to introduce the saturation (one that takes into account the relative inclinations of the field and the flow, e.g. the cross-convergence). Lastly, more realistic boundary conditions, and a parameter study with increasing Rem will be performed. The final goal is to be able to identify the more of operation during the (non-linear) saturation phase of the dynamo at high Rem in order to predict the likely topology of the magnetic field that one might observe at stars such as Betelgeuse.

#### Acknowledgments.

SBFD thanks the LOC of IAU Symposium No. 210 for support. The radiation-hydrodynamic simulations yielding the input flow field used in this study was computed at the computing center of the university in Kiel, at the Theoretical Astrophysics Center in Copenhagen, and on the Sun cluster at Ångström laboratory in Uppsala. The dynamo calculations were performed at the Institute for Solar Physics of the Royal Swedish Academy of Sciences.

## References

Archontis, V.D., Dorch, S.B.F., & Nordlund, Å. 2002, A&A, accepted

Archontis, V.D., & Dorch, S.B.F. 1999, in Stellar Dynamos: Non-linearity and chaotic flows, ASP conf. series 178, 1

Brandenburg, A. 2001, ApJ, 550, 824

Cattaneo, F. 1999, ApJ, 515, L39

Dorch, S.B.F. 2000, Physica Scripta, 61, 717

Dorch, S.B.F., & Ludwig, H.-G. 2002, Astron. Nachr., 323, in press

Dorch, S.B.F., & Nordlund, Å. 2001, A&A, 365, 562

Freytag, B., & Finnsson, S. 2002, ibid

Freytag, B., & Mizuno-Wiedner, M. 2002, ibid

Freytag, B., Steffen, M., & Dorch, S.B.F. 2002, Astron. Nachr., 323, 213

Galsgaard, K., & Nordlund, Å. 1997, Journ. Geophysical Res., 101, 219

Schrijver, C., & Zwaan, C. 2000, Solar and Stellar Magnetic Activity, CAPS, 34, Cambridge

Stein, R.F. 2002, ibid