format: text

text:
# 𝜌→𝜋𝜋 decay in nuclear medium

## 1 Introduction

In recent years a lot of efforts have been undertaken in order to understand the properties of hadrons in nuclear medium. This challenging theoretical issue gains a lot of importance at the beginning of the operation of RHIC, where a proper and accurate inclusion of hadrons is necessary to describe the evolution of the “hadronic soup” formed in the collision. It is commonly accepted that hadrons are significantly modified by nuclear matter. The arguments range from simple scaling of masses , through numerous hadronic model calculations , QCD sum-rule techniques , approaches motivated by the chiral symmetry , to model-independent predictions based on low-density expansion and dispersion relations . In these calculations masses of hadrons, or widths, or both, are significantly changed by the presence of the medium (for a recent review see ). In addition, medium may induce meson mixing absent in the vacuum . An indirect indication for the modification of vector mesons has been provided by the dilepton production measurements in relativistic heavy-ion collisions (CERES , HELIOS ). The dilepton productions has been theoretically studied with medium-modified vector mesons , which helps to explain the low-mass enhancement of the dilepton yields.

There is a lot of studies of the meson two-point functions in the literature, but only a few devoted to meson three-point functions. In their study of the $`\rho `$-meson in-medium spectral function Herrmann, Friman, and Nörenberg  have computed the $`\rho \pi \pi `$ vertex for the $`\rho `$ at rest in nuclear matter. They have applied a hadronic model with the $`\mathrm{\Delta }`$ isobar and with non-relativistic couplings. Temperature effects on the $`\rho \pi \pi `$ interaction have been considered by Song and Koch . Urban, Buballa, Rapp and Wambach  have extended the calculation of Ref.  to non-zero three-momenta of the $`\rho `$ meson and non-zero temperature. The authors of the present paper have analyzed the $`\omega \pi \pi `$ decay in nuclear medium . Krippa has theoretically studied the effects of density on chiral mixing of meson three-point functions . In Ref.  the $`\rho \pi \pi `$ coupling has been studied in quark matter. Otherwise, the topic of medium effects on hadronic couplings is very much terra incognita. Certainly, if the two-point functions are significantly altered by the medium, one expects that the three-point functions should also change. The issue is important for modeling the hadronic evolution in relativistic heavy-ion collisions, since the change of the hadron couplings results in altered transition rates between hadrons.

Out of many meson couplings, the $`\rho \pi \pi `$ vertex is especially important, since the $`\rho `$ plays an essential role in hadron dynamics. The vacuum value of the coupling constant is large, $`g_{\rho \pi \pi }6`$. Since the couplings of $`\rho `$ and pions to nucleons and $`\mathrm{\Delta }`$ isobars are also large, we expect significant medium modifications of $`g_{\rho \pi \pi }`$. As our calculation shows, this is indeed the case, with the coupling significantly altered already at the nuclear saturation density. The effect is mainly due to the interactions with the $`\mathrm{\Delta }`$. Our method is similar to Refs. , with the following differences: We use a fully relativistic framework, with relativistic interactions. The leading-density approximation is used, which makes the calculation very simple, as no integration over the nucleon momenta is necessary. We analyze analyticity in the virtual $`\rho `$ mass, which is non-trivial in our problem due to low-mass thresholds. These thresholds largely influence the results. In addition, keeping non-zero three-momentum of the $`\rho `$, as in Ref. , allows us to look separately on the longitudinal and transverse polarizations. In our comparison to the CERES dilepton data we use the fire-cylinder expansion model of Refs.  and take into account the experimental cuts, which are very important for the detailed numerical analysis.

The outline and the main results of the paper are as follows: In Sect. 2 we introduce our framework, i.e. the relativistic model of mesons interacting with nucleons and $`\mathrm{\Delta }`$ isobars, the latter treated as Rarita-Schwinger fields. We introduce the necessary hadronic vertices and set the coupling constants. As already mentioned, we work at zero temperature and to the leading-order in baryon density. In Sect. 3 we present our results for the $`\rho \pi \pi `$ vertex for the $`\rho `$ meson at rest with respect to the medium, and point out large medium effects even at the nuclear saturation density. We investigate analyticity of the coupling in the virtual $`\rho `$ meson mass and point out its significance for the quantitative results. We look at effects of the energy-dependent width of the $`\mathrm{\Delta }`$ and find they are not very large. In Sect. 4 we analyze the general case of finite three-momentum of the $`\rho `$ with respect to the medium, and calculate the widths and spectral functions in the transverse and longitudinal channels. In Sect. 5 we apply the obtained spectral functions to compute the dilepton production rate from $`\rho `$ decays for the CERES $`Pb+Au`$ experiment. The calculation carefully includes the CERES experimental kinematic cuts, which is crucial for the results. We also include the effects of the expansion of the fire cylinder, which enter at the level of a few percent. It is found that spreading of the width helps to explain the low-mass enhancement of the dilepton spectra, but the overall normalization fails significantly short of the experimental data. Section 6 presents our conclusions and a discussion of several additional points.

## 2 The framework

Throughout this paper we choose conventionally $`q`$ as the incoming momentum of the $`\rho `$ carrying isospin index $`b`$ and polarization vector $`\epsilon ^\mu `$, $`p`$ as the outgoing momentum of the pion carrying isospin index $`a`$, and $`qp`$ as the outgoing momentum of the pion carrying isospin index $`c`$. With this convention, the vacuum value of the $`\rho \pi \pi `$ vertex (Feynman rule) is

$$iV_{\rho _\mu ^b\pi ^a\pi ^c}=g_\rho ϵ^{acb}(2p^\mu q^\mu ).$$
(1)
The medium modifications of the coupling, calculated below, will be compared to Eq. (1).

Our calculation of the in-medium $`\rho \pi \pi `$ vertex is made in the framework of a fully relativistic hadronic theory, where mesons interact with the nucleons and $`\mathrm{\Delta }`$ isobars. Other in-medium calculations indicate that the leading-density approximation is sufficient up to densities of the order of the nuclear saturation density.333An extension beyond the leading-density approximation would require accounting for the Fermi-sea motion, which technically leads to keeping the integration over $`𝐤`$ in our expressions, and, more importantly, the inclusion of correlations. In many similar calculations the Fermi-motion effects have been found to influence the results weakly at moderate densities. The inclusion of nucleon correlations is a difficult problem, extending far beyond the present work. Also, in Ref.  the finite-temperature effects in the $`\rho \pi \pi `$ coupling have been found to be moderate for temperatures up to about $`150`$MeV, typical for relativistic heavy-ion collisions (cf. Fig. 7(a,b) of ). This justifies the use of zero temperature to obtain the vertices of Fig. 1, at least as a first approximation.

To the leading-density order only the diagrams shown in Fig. 1 contribute to the $`\rho \pi \pi `$ process. The solid lines denote $`i`$ times the in-medium nucleon propagator, which can be decomposed in the usual way into the free and density parts :

$`iG(k)`$ $``$ $`iG_F(k)+iG_D(k)=`$
$`=`$ $`i(/k+m_N)\left[{\displaystyle \frac{1}{k^2m_N^2+i\epsilon }}+{\displaystyle \frac{i\pi }{E_k}}\delta (k_0E_k)\theta (k_F|𝐤|)\right],`$
where $`k`$ denotes the nucleon four-momentum, $`m_N`$ is the nucleon mass, $`E_k=\sqrt{m_N^2+𝐤^2}`$, and $`k_F`$ is the Fermi momentum of nuclear matter. The double line in the diagrams of Fig. 1 denotes $`i`$ times the relativistic $`\mathrm{\Delta }`$ propagator,

$`iG_\mathrm{\Delta }^{\alpha \beta }(k)`$ $`=`$ $`i{\displaystyle \frac{/k+m_\mathrm{\Delta }}{k^2\left(m_\mathrm{\Delta }\frac{i}{2}\mathrm{\Gamma }_\mathrm{\Delta }\right)^2}}\left(g^{\alpha \beta }+{\displaystyle \frac{1}{3}}\gamma ^\alpha \gamma ^\beta +{\displaystyle \frac{2k^\alpha k^\beta }{3m_\mathrm{\Delta }^2}}+{\displaystyle \frac{\gamma ^\alpha k^\beta \gamma ^\beta k^\alpha }{3m_\mathrm{\Delta }}}\right).`$
This formula corresponds to the usual Rarita-Schwinger definition , with the parameter choice $`A=1`$. We have modified the denominator in Eq. (LABEL:Dprop) in order to account for the finite width of the $`\mathrm{\Delta }`$ resonance.

Since we are interested in density effects, one of the nucleon lines in each of the diagrams of Fig. 1 must involve the nucleon density propagator, $`G_D`$. For kinematic reasons, diagrams with more than one $`G_D`$ vanish. The wavy line in Fig. 1 denotes the $`\rho `$ meson, and the dashed lines correspond to pions. All external particles are on mass shell.

The meson-nucleon vertices (Feynman rules) needed for our calculation have the standard form

$`iV_{\pi ^aNN}`$ $`=`$ $`{\displaystyle \frac{g_A}{2F_\pi }}/p\gamma _5\tau ^a,`$ (4)
$`iV_{\rho _\mu ^bNN}`$ $`=`$ $`ig_\rho (\gamma ^\mu +{\displaystyle \frac{i\kappa _\rho }{2m_N}}\sigma ^{\mu \nu }q_\nu ){\displaystyle \frac{\tau ^b}{2}},`$ (5)
$`iV_{\rho _\mu ^b\pi ^aNN}`$ $`=`$ $`i{\displaystyle \frac{g_\rho g_A}{2F_\pi }}\gamma ^\mu \gamma _5\epsilon ^{abc}\tau _c`$ (6)
where $`p`$ is the outgoing four-momentum of the pion, $`q`$ is the incoming momentum of the $`\rho `$ meson, and $`a`$ and $`b`$ are the isospin indices of the pion and the $`\rho `$, respectively. We have chosen the pseudovector pion-nucleon coupling. The vertex (6) follows from the minimal substitution in Eq. (4).

Whereas the couplings of mesons to nucleons listed above are well established and the corresponding parameters are well known (except for the long-lasting controversy with $`\kappa _\rho `$), the relativistic couplings involving the $`\mathrm{\Delta }`$ resonance are a topic of an on-going research and discussion . Several structures are possible for these vertices, and ambiguities related to the choice of the so-called off-shell parameters have not been resolved. For our pragmatic goal of estimating the size of the in-medium $`\rho \pi \pi `$ vertex, we will adopt the following simple-minded and popular philosophy: all off-shell couplings are set to zero.

The meson-nucleon-$`\mathrm{\Delta }`$ vertices have the following form

$`iV_{\pi ^aN\mathrm{\Delta }_\alpha }`$ $`=`$ $`{\displaystyle \frac{f_{\pi N\mathrm{\Delta }}}{m_\pi }}\theta ^{\alpha \nu }(Z)p_\nu T^a,`$ (7)
$`iV_{\rho _\mu ^bN\mathrm{\Delta }_\alpha }`$ $`=`$ $`ig_1(\theta ^{\mu \alpha }(Y)/q\gamma _5\theta ^{\mu \nu }(Y)q_\nu \gamma ^\alpha \gamma _5)T^b`$ (8)
$`+ig_2(\theta ^{\mu \alpha }(X)\gamma _5kq\theta ^{\mu \nu }(X)q_\nu \gamma _5k^\alpha )T^b,`$
$`iV_{\rho _\mu ^b\pi ^aN\mathrm{\Delta }_\alpha }`$ $`=`$ $`i{\displaystyle \frac{g_\rho f_{\pi N\mathrm{\Delta }}}{m_\pi }}\theta ^{\alpha \mu }(Z)\epsilon ^{abc}T_c,`$ (9)
where $`k`$ is the nucleon four momentum, and

$$\theta ^{\sigma \lambda }(W)=g^{\sigma \lambda }(W+\frac{1}{2})\gamma ^\sigma \gamma ^\lambda ,W=Z,Y,X.$$
(10)
The combination $`(W+\frac{1}{2})`$ is called the off-shell parameter. According to the prescription stated above, we set $`Z=Y=X=\frac{1}{2}`$, such that all off-shell parameters vanish, and we simply have $`\theta ^{\sigma \lambda }(W)=g^{\sigma \lambda }`$. Furthermore, we take arbitrarily $`g_2=0`$. The matrices $`T^a`$ in Eqs. (7-9) are the standard isospin $`\frac{1}{2}\frac{3}{2}`$ transition matrices, given in App. A. The vertex (9) follows from the minimum substitution in Eq. (7).

For the $`\pi \mathrm{\Delta }\mathrm{\Delta }`$ coupling there are, according to Ref. , three possible structures:

$$iV_{\pi ^a\mathrm{\Delta }_\alpha \mathrm{\Delta }_\beta }=\left(G_1g^{\alpha \beta }/p\gamma _5+G_2(\gamma ^\alpha p^\beta +p^\alpha \gamma ^\beta )\gamma _5+G_3\gamma ^\alpha /p\gamma _5\gamma ^\beta \right)T_\mathrm{\Delta }^a.$$
(11)
We drop the off-shell couplings by setting $`G_2=G_3=0`$. For the $`\rho \mathrm{\Delta }\mathrm{\Delta }`$ vertex we use the minimal vector current coupling  and the universality prescription, which gives

$$iV_{\rho _\mu ^b\mathrm{\Delta }_\alpha \mathrm{\Delta }_\beta }=ig_\rho (\gamma ^\mu g^{\alpha \beta }+g^{\alpha \mu }\gamma ^\beta +g^{\beta \mu }\gamma ^\alpha +\gamma ^\alpha \gamma ^\mu \gamma ^\beta )T_\mathrm{\Delta }^b.$$
(12)
The constants $`g_1`$ and $`G_1`$ are adjusted in such a way, that in the nonrelativistic limit we recover the couplings $`\sqrt{2}(f^{}/m_\pi )\epsilon _{ijk}p^iS_\mathrm{\Delta }^jT_\mathrm{\Delta }^a`$ and $`(f_\mathrm{\Delta }/m_\pi )p_jS_\mathrm{\Delta }^jT_\mathrm{\Delta }^a`$, respectively , where $`k`$ is the spin index of the $`\rho `$. The comparison, with the explicit form of the Rarita-Schwinger spinors and the matrices $`T_\mathrm{\Delta }`$ and $`S_\mathrm{\Delta }`$ (see App. B) gives $`G_1=`$ $`\frac{3}{2}f_\mathrm{\Delta }/m_\pi `$ and $`g_1=`$ $`\sqrt{2}f^{}/m_\pi `$. Our choice of the physical parameters is as follows :

$`g_A`$ $`=`$ $`1.26,F_\pi =93\mathrm{M}\mathrm{e}\mathrm{V},m_\pi =139.6\mathrm{MeV},`$
$`g_\rho `$ $`=`$ $`5.26,\kappa _\rho =6,`$ (13)
$`f_{\pi N\mathrm{\Delta }}`$ $`=`$ $`2.12,f^{}=2.12,f_\mathrm{\Delta }=0.802.`$
In addition to the mentioned ambiguities in the choice of the form of the $`\mathrm{\Delta }`$ couplings , there are uncertainties in the values of the coupling constants. For some constants one typically uses the quark model predictions, or the large-$`N_c`$ arguments which relate the $`\mathrm{\Delta }`$ couplings to the nucleon couplings . Any adopted scheme should fit the values of the coupling constants to the available data for various processes. However, for our purpose of estimating the size of the medium effect on the $`\rho \pi \pi `$ vertex these ambiguities are not essential. In addition, the result, as we shall see shortly, is dominated by diagram (g) of Fig. 1, containing only the $`\pi N\mathrm{\Delta }`$ and $`\pi \rho N\mathrm{\Delta }`$ couplings which are well established. Our results are not very sensitive to the choice of other couplings. Also, for simplicity of the approach and from the lack of knowledge we do not include any form-factors in the vertices.

We work in the rest frame of the nuclear matter, however effort is made to write all expressions covariantly, which turns out to be very convenient. Our calculation is made in the following way. First, we evaluate the diagrams of Fig. 1. The result for the full $`\rho \pi \pi `$ amplitude has the generic form $`_{acb}=\epsilon _\mu A_{acb}^\mu `$, where the vertex function is

$`A_{acb}^\mu `$ $`=`$ $`ϵ^{acb}(A_{\mathrm{vac}}^\mu +A_{\mathrm{med}}^\mu ),`$ (14)
$`A_{\mathrm{vac}}^\mu `$ $`=`$ $`g_\rho (2p^\mu q^\mu ),`$
$`A_{\mathrm{med}}^\mu `$ $`=`$ $`{\displaystyle \frac{d^3k}{(2\pi )^3}\frac{m_N}{E_k}(Ap^\mu +Bq^\mu +Ck^\mu )\theta (k_F|𝐤|)},`$
and $`A`$, $`B`$, and $`C`$ are scalar functions depending on scalar products of the four-vectors $`q`$, $`p`$, and $`k`$, with $`k^0=E_k`$. The occupation function is made explicitly Lorentz-invariant when we write $`|𝐤|=\sqrt{(ku)^2m_N^2}`$, where $`u`$ is the four-velocity of the medium. The term with $`k^\mu `$, upon the evaluation of the integral, can be in general proportional to the three Lorentz vectors present in the problem, namely

$$\frac{d^3k}{(2\pi )^3}\frac{m_N}{E_k}Ck^\mu \theta (k_F|𝐤|)=C_pp^\mu +C_qq^\mu +C_uu^\mu ,$$
(15)
where $`C_q`$, $`C_p`$, and $`C_u`$ are scalar functions of $`p^2`$, $`q^2`$, $`pq`$, $`pu`$, $`qu`$, and $`k_F`$. Contracting Eq. (15) with $`p_\mu `$, $`q_\mu `$, and $`u_\mu `$ we obtain a set of linear algebraic equations for $`C_q`$, $`C_p`$, and $`C_u`$, which can be solved. However, at the leading-density approximation the problem becomes even simpler. We can work in the rest frame of the medium, where $`u^\mu =(1,0,0,0)`$. It is obvious that the leading-density approximation is equivalent to setting the three-vector $`𝐤`$ to zero in the functions $`A`$, $`B`$, $`C`$, and $`E_k`$ appearing in the integrands of Eqs. (14,15). Then $`\frac{d^3k}{(2\pi )^3}\theta (k_f|𝐤|)=\frac{1}{4}\rho _B`$, while higher-order terms in $`𝐤`$ generate terms with higher exponents of the baryon density $`\rho _B`$. Now, with $`𝐤=0`$, $`k^0=m_N`$, the contraction of Eq. (15) with $`q_\mu `$, $`p_\mu `$, and $`u_\mu `$ gives the set of equations

$`{\displaystyle \frac{1}{4}}\rho _Bm_Np^0\overline{C}`$ $`=`$ $`C_pp^2+C_qqp+C_up^0,`$
$`{\displaystyle \frac{1}{4}}\rho _Bm_Nq^0\overline{C}`$ $`=`$ $`C_ppq+C_qq^2+C_uq^0,`$ (16)
$`{\displaystyle \frac{1}{4}}\rho _Bm_N\overline{C}`$ $`=`$ $`C_ppu+C_qqu+C_u,`$
where $`\overline{C}`$ is $`C`$ with $`𝐤=0`$. Since in the general case vectors $`p`$, $`q`$, and $`u`$ are linearly-independent, the solution of Eqs. (16) is $`C_p=C_q=0`$, $`C_u=\frac{1}{4}\rho _Bm_N\overline{C}`$. Thus, only the term proportional to $`u^\mu `$ in Eq. (15) is present in the leading-density approximation. In other words, we have

$$A_{\mathrm{med}}^\mu =\frac{1}{4}\rho _B(\overline{A}p^\mu +\overline{B}q^\mu +\overline{C}m_Nu^\mu ),$$
(17)
where the coefficients $`\overline{A}`$, $`\overline{B}`$, and $`\overline{C}`$ are obtained from $`A`$, $`B`$, and $`C`$ by simply setting $`𝐤=0`$.

To summarize this part, we restate the necessary steps needed to obtain the leading-density amplitude: The traces in diagrams of Fig. 1 are evaluated, leading to Eq. (14). Then we replace $`k^\mu `$ by $`m_Nu^\mu `$, set $`𝐤=0`$, and arrive at Eq. (17). Certainly, this very simple method is general for any problem involving baryon loops with density-dependent nucleon propagators. In our calculation we have used a standard Dirac algebra package . The isospin traces are evaluated in App. B.

## 3 Results for $`\rho `$ decaying at rest

We begin the presentation of the results with the case where the $`\rho `$ is at rest with respect to nuclear matter, $`𝐪=0.`$ In this kinematics we find $`\overline{B}=\frac{1}{2}\overline{A}`$, and $`\overline{C}=0`$. The fact that $`\overline{C}=0`$ is reflecting the equality of the transversely and longitudinally polarized $`\rho `$ meson propagators at $`𝐪=0`$, as will be shown in Sect. 4. The result $`\overline{B}=\frac{1}{2}\overline{A}`$ is consistent with the Ward-Takahashi identity $`q_\mu A^\mu =D_\pi ^1(p)D_\pi ^1(pq)`$, where $`D_\pi `$ denotes the pion propagator dressed with nucleon and nucleon-$`\mathrm{\Delta }`$ bubbles . With $`𝐪=0`$ and pions on the mass shell, we obviously have $`D_\pi ^1(p)=D_\pi ^1(pq)`$, and, consequently, $`\overline{A}pq+\overline{B}q^2=0`$, which immediately gives $`\overline{B}=\frac{1}{2}\overline{A}`$. Therefore, for $`𝐪=0`$ the in-medium vertex function is proportional to $`2p^\mu q^\mu `$:

$$A_{\mathrm{med}}^\mu (𝐪=0)=\frac{1}{8}\rho _B\overline{A}(𝐪=0)(2p^\mu q^\mu ).$$
(18)
The full formula for $`\overline{A}(𝐪=0)`$ is very lengthy, hence we present below only the contribution from the bubble diagrams (f,g), which are simple. In App. C we list the contributions to $`\overline{A}(𝐪=0)`$ from all diagrams in the special case of $`m_\pi =0`$.444It should be stressed that there is no obvious expansion parameter in the problem. Expanding the amplitude in $`m_\pi `$ or $`m_\mathrm{\Delta }m_N`$ and keeping the lowest terms does not lead to a good approximation. It turns out that the $`N\mathrm{\Delta }`$ bubble diagram (g) is the dominant one. We have

$`\overline{A}^{(g)}(𝐪=0)`$ $`=`$ $`{\displaystyle \frac{16}{9}}g_\rho \left({\displaystyle \frac{f_{\pi N\mathrm{\Delta }}}{m_\pi }}\right)^2\times `$
$`{\displaystyle \frac{2(m_\mathrm{\Delta }^2m_\pi ^2\frac{1}{2}m_NM)(m_N+m_\mathrm{\Delta }+\frac{1}{2}M)}{m_\mathrm{\Delta }^2(m_N^2(m_\mathrm{\Delta }\frac{i}{2}\mathrm{\Gamma }_\mathrm{\Delta })^2+m_\pi ^2+m_NM)}}+(MM),`$
with $`M`$ denoting the mass of the $`\rho `$ meson. For comparison, the contribution from the $`NN`$ bubble, Fig. 1(f), is

$$\overline{A}^{(f)}(𝐪=0)=g_\rho \left(\frac{g_A}{2F_\pi }\right)^2\frac{32m_Nm_\pi ^2}{m_\pi ^4m_N^2M^2}.$$
(20)
In the limit of large $`m_N`$, with $`m_\mathrm{\Delta }m_N`$ fixed, and $`\mathrm{\Gamma }_\mathrm{\Delta }=0`$, expression (3) reduces to

$$\overline{A}^{(g)}(𝐪=0)\frac{16}{9}g_\rho \left(\frac{f_{\pi N\mathrm{\Delta }}}{m_\pi }\right)^2\frac{4(m_\mathrm{\Delta }m_N)}{M^2/4(m_\mathrm{\Delta }m_N)^2},$$
(21)
which agrees with non-relativistic calculations.

In the following we shall treat $`M`$ as the mass of a virtual $`\rho `$ meson. Virtual $`\rho `$ mesons are needed for the analysis of the dilepton production in Sect. 5. Analyticity of the vertex function is nontrivial in the variable $`M`$. We can see from the denominators of Eqs. (3,20), that for $`\mathrm{\Gamma }_\mathrm{\Delta }=0`$ singularities occur at

$`M^2`$ $`=`$ $`\left({\displaystyle \frac{m_\mathrm{\Delta }^2m_N^2m_\pi ^2}{m_N}}\right)^2=(0.657\mathrm{GeV})^2,`$ (22)
$`M^2`$ $`=`$ $`\left({\displaystyle \frac{m_\pi ^2}{m_N}}\right)^2=(0.021\mathrm{GeV})^2.`$ (23)
Triangle diagrams of Fig. 1 also have the above singularities, and additionally bring in high-lying particle-antiparticle production singularities at $`M^2=(2m_N)^2`$ and at $`M^2=\left(m_\mathrm{\Delta }+m_N\right)^2`$, which are physically not relevant. The above analytic structure is manifest in the numerical calculations presented below. For non-zero $`\mathrm{\Gamma }_\mathrm{\Delta }`$ the pole at (22) changes to a broader structure. Thus, analyticity is important — it immediately leads to large changes of the vertex function near the poles.

Since for the case of $`𝐪=0`$ the matter-induced vertex function (18) has the same Lorentz structure as in the vacuum, i.e. proportional to $`2p^\mu q^\mu `$, it is convenient for our quantitative studies to introduce an effective $`\rho \pi \pi `$ coupling constant, defined as

$$g_{\mathrm{eff}}=\left|g_\rho +\frac{1}{8}\rho _B\overline{A}(𝐪=0)\right|.$$
(24)
The absolute value is taken, since with non-zero $`\mathrm{\Gamma }_\mathrm{\Delta }`$ the quantity $`\overline{A}(𝐪=0)`$ is complex. In Fig. 2 we plot the ratio $`g_{\mathrm{eff}}/g_\rho `$ as a function of the virtual $`\rho `$ mass, $`M`$, at the saturation density, $`\rho _B=\rho _0`$. The dashed line corresponds to the zero-width $`\mathrm{\Delta }`$. We can clearly see the singularities of Eq. (22,23), whose positions are indicated by arrows. The solid line in Fig. 2 shows the calculation with the vacuum value of the $`\mathrm{\Delta }`$ width, $`\mathrm{\Gamma }_\mathrm{\Delta }=120\mathrm{M}\mathrm{e}\mathrm{V}.`$ In this case the pole at (22) changes to a broad structure. The pole of Eq. (23) remains, of course, unchanged. The considerable difference between the solid and dashed curves in the range of $`M`$ between $`0.6`$ and $`1`$GeV shows that the results are sensitive to the assumed value for $`\mathrm{\Gamma }_\mathrm{\Delta }.`$ We note that at low masses $`M`$, between $`0.07`$ and $`0.55\mathrm{GeV}`$, the effective coupling $`g_{\mathrm{eff}}`$ is lower than the vacuum value, thus the medium lowers the coupling, while above $`M0.55\mathrm{GeV}`$ the effect is opposite: the coupling is increased. Around the physical $`\rho `$ mass, $`M=m_\rho `$, the effective coupling is roughly two times larger than in the vacuum. For the width of the $`\rho \pi \pi `$ decay this means a factor of $`4`$ enhancement, giving an in-medium width to the $`\rho `$ of about $`600`$MeV at the saturation density. A similar estimate has been obtained e.g. in Refs. .

The dot-dashed line in Fig. 2 shows the result of a calculation with the energy-dependent width, parameterized as a function of the $`s`$ variable as follows:

$`\mathrm{\Gamma }_\mathrm{\Delta }(s)`$ $`=`$ $`\mathrm{\Gamma }_\mathrm{\Delta }\left({\displaystyle \frac{q_{\mathrm{cm}}(s)}{q_{\mathrm{cm}}(m_\mathrm{\Delta })}}\right)^3\theta (s(m_N+m_\pi )^2),`$ (25)
$`q_{\mathrm{cm}}(s)`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{s}}}\sqrt{(s(m_N+m_\pi )^2)(s(m_Nm_\pi )^2)}.`$ (26)
As can be seen from Fig. 2, the effects of energy-dependent widths are not large in our analysis.

It is interesting to look at the anatomy of the relative medium contribution to $`g_{\mathrm{eff}}`$, i.e.  of the quantity $`\frac{1}{8}\rho _B\overline{A}(𝐪=0)/g_\rho `$. At $`M=m_\rho =776`$MeV and $`\rho _B=\rho _0=0.17\mathrm{fm}^3`$ we find that the diagrams (a)-(g) contribute, correspondingly, $`0.08`$, $`0.220.09i`$, $`0.0080.008i`$, $`0.100.11i`$, $`0.11+0.63i`$, $`0.008`$, and $`0.621.16i`$, with the total of $`0.730.75i`$, for the case $`\mathrm{\Gamma }_\mathrm{\Delta }=120\mathrm{M}\mathrm{e}\mathrm{V}`$, and $`.08`$, $`0.33`$, $`0.001`$, $`0.25`$, $`0.87`$, $`0.008`$, and $`2.09`$, with the total of $`1.71`$, for the case $`\mathrm{\Gamma }_\mathrm{\Delta }=0`$. As advocated above, the largest contribution comes from the $`N\mathrm{\Delta }`$ bubble diagram (g).

In Fig. 3 we display the results for different values of the baryon density, $`\rho _B`$. In this study we have neglected the energy-dependence of the $`\mathrm{\Delta }`$ width, but included the density dependence of $`\mathrm{\Gamma }_\mathrm{\Delta }`$. The $`\mathrm{\Delta }`$ broadens moderately at the nuclear saturation density , with Pauli blocking giving a decrease of $`\mathrm{\Gamma }_\mathrm{\Delta }`$ by about 40MeV, and absorption processes giving an increase by about 80MeV, such that the net effect is an increase by about 40MeV. Therefore we parameterize

$$\mathrm{\Gamma }_\mathrm{\Delta }(\rho _B)=\mathrm{\Gamma }_\mathrm{\Delta }(1+\frac{40\mathrm{M}\mathrm{e}\mathrm{V}}{\mathrm{\Gamma }_\mathrm{\Delta }}\frac{\rho _B}{\rho _0}).$$
(27)
Obviously, the effective coupling increases with $`\rho _B`$, as can be seen in Fig. 3. A caveat is in place here. Our method can be trusted numerically only at low values of the baryon density, such that the leading-density approximation holds. On the other hand we can see from Fig. 3 that the effects are large already at the saturation density, and certainly the approximation breaks down at larger values of $`\rho _B`$. Therefore our numerical results at high densities, here and in the following parts of the paper, have to be taken with a grain of salt and treated as indication of possible large effects rather than accurate numerical predictions.

Finally, we complete our discussion of the in-medium $`\rho \pi \pi `$ vertex at $`𝐪=0`$ with Fig. 4, which shows the result of the calculation with fixed $`\mathrm{\Gamma }_\mathrm{\Delta }`$, but with $`m_N`$ and $`m_\mathrm{\Delta }`$ scaled down to 70% of their vacuum values. A decrease of that order at the saturation density is anticipated from several approaches . We notice that $`g_{\mathrm{eff}}`$ is enhanced and shifted to lower values of $`M`$ when the baryon masses are rescaled.

## 4 Results for moving $`\rho `$

When the $`\rho `$ meson moves with a non-zero momentum $`𝐪`$ with respect to the medium, its propagation is different for transverse and longitudinal polarizations, defined by quantizing the spin along the direction of $`𝐪`$. To analyze this effect we shall consider the width of the transversely and longitudinally polarized $`\rho `$ mesons due to the decay into two pions. This allows us to present the result in a more compact form, rather than looking separately at the functions $`\overline{A}`$, $`\overline{B}`$, and $`\overline{C}`$. The expression for the width for the decay $`\rho ^0\pi ^+\pi ^{}`$, as viewed from the rest frame of the medium, is

$$\mathrm{\Gamma }_{\rho \pi \pi }=\frac{1}{n_s}\underset{s}{}\frac{1}{2q_0}\frac{d^3p}{(2\pi )^32p_0}\frac{d^3p^{}}{(2\pi )^32p_0^{}}||^2(2\pi )^4\delta ^{(4)}(qpp^{}),$$
(28)
where $`n_s`$ is the number of spin states of the $`\rho `$ meson, and $`_s`$ denotes the sum over these spin states. The division by $`q_0=\sqrt{M^2+𝐪^2}`$ in Eq. (28), rather than by $`M`$, accounts for the time-dilatation effect. We perform the phase-space integral in the rest frame of the nuclear medium, and obtain

$$\mathrm{\Gamma }_{\rho \pi \pi }=\frac{1}{n_s}\underset{s}{}\frac{1}{2q_0}\underset{b=1,2}{}\underset{0}{\overset{\gamma ^{}}{}}\mathrm{sin}\gamma \frac{(𝐩^{(b)})^2}{8\pi p_0^{(b)}(q_0p_0^{(b)})|a^{(b)}|}||^2d\gamma ,$$
(29)
where $`_b`$ is the sum over the two possible kinematic branches, which can appear when the two-body decay is viewed from a frame where the decaying particle moves. The angle $`\gamma `$ is between the directions of $`𝐪`$ and $`𝐩`$. The second branch appears only for $`|𝐪|`$ above a critical value,

$$|𝐪|>\frac{M\sqrt{M^24m_\pi ^2}}{2m_\pi }q_{\mathrm{crit}}.$$
(30)
Elementary kinematic considerations give

$`|𝐩^{(1,2)}|`$ $`=`$ $`{\displaystyle \frac{M^2|𝐪|\mathrm{cos}\gamma \pm q_0\sqrt{M^44m_\pi ^2(M^2+𝐪^2\mathrm{sin}^2\gamma )}}{2(M^2+𝐪^2\mathrm{sin}^2\gamma )}},`$
$`q_0`$ $`=`$ $`\sqrt{M^2+𝐪^2},p_0^{(1,2)}=\sqrt{m_\pi ^2+\left(𝐩^{(1,2)}\right)^2},`$
$`a^{(1,2)}`$ $`=`$ $`{\displaystyle \frac{d(q_0\sqrt{m_\pi ^2+r^2}\sqrt{m_\pi ^2+r^22r|𝐪|\mathrm{cos}\gamma +𝐪^2})}{dr}}|_{r=|𝐩^{(1,2)}|},`$
$`\gamma ^{}`$ $`=`$ $`\{\begin{array}{c}\pi \mathrm{for}qq_{\mathrm{crit}}\hfill \\ \mathrm{arcsin}\left(\frac{M\sqrt{M^24m_\pi ^2}}{2m_\pi |𝐪|}\right)\mathrm{for}q>q_{\mathrm{crit}}\hfill \end{array}.`$ (33)
The transversely polarized $`\rho `$ has two helicity states ($`n_s=2`$), with projection $`s=\pm 1`$ on the direction of $`𝐪,`$ described by polarization vectors $`\epsilon _{(\pm )}^\mu `$, while the longitudinally polarized $`\rho `$ has one helicity state ($`n_s=1`$), with the corresponding projection $`s=0`$, described by the polarization vector $`\epsilon _{(0)}^\mu `$. An explicit calculation yields

$`\epsilon _{(+)}^\mu \epsilon _{(+)}^\nu \epsilon _{()}^\mu \epsilon _{()}^\nu `$ $`=`$ $`g^{\mu \nu }u^\mu u^\nu {\displaystyle \frac{(q^\mu quu^\mu )(q^\nu quu^\nu )}{qq(qu)^2}}T^{\mu \nu },`$
$`\epsilon _{(0)}^\mu \epsilon _{(0)}^\nu `$ $`=`$ $`{\displaystyle \frac{q^\mu q^\nu }{qq}}+u^\mu u^\nu +{\displaystyle \frac{(q^\mu quu^\mu )(q^\nu quu^\nu )}{qq(qu)^2}}L^{\mu \nu }.`$
Note that by summing over all polarizations one recovers the usual formula, i.e.

$$T^{\mu \nu }+L^{\mu \nu }=g^{\mu \nu }\frac{q^\mu q^\nu }{qq}.$$
(35)
The tensors $`T^{\mu \nu }`$ and $`L^{\mu \nu }`$ are defined with such signs as to form projection operators, i.e. , $`T^{\mu \nu }T_\nu ^\alpha =T^{\mu \alpha }`$, $`L^{\mu \nu }L_\nu ^\alpha =L^{\mu \alpha }`$, and $`T^{\mu \nu }L_\nu ^\alpha =0`$. Furthermore, we have $`T^{\mu \nu }q_\nu =0`$ and $`L^{\mu \nu }q_\nu =0`$, which reflects current conservation, as well as $`T^{\mu \nu }u_\nu =0`$. Through the use of these relation and Eqs. (14,17) we find that

$`|_T|^2={\displaystyle \underset{s=\pm }{}}\epsilon _{(s)}^\mu A_\mu ^{}\epsilon _{(s)}^\nu A_\nu =\left|2(g_\rho +{\displaystyle \frac{1}{8}}\rho _B\overline{A})\right|^2p_\mu T^{\mu \nu }p_\nu ,`$ (37)
$`|_L|^2={\displaystyle \underset{s=0}{}}\epsilon _{(s)}^\mu A_\mu ^{}\epsilon _{(s)}^\nu A_\nu =`$
$``$ $`(2(g_\rho +{\displaystyle \frac{1}{8}}\rho _B\overline{A})^{}p_\mu +\overline{C}^{}m_Nu_\mu )L^{\mu \nu }(2(g_\rho +{\displaystyle \frac{1}{8}}\rho _B\overline{A})p_\nu +\overline{C}m_Nu_\nu ).`$
Note that the value of the coefficient $`\overline{B}`$ in Eq. (17) is irrelevant for the widths. Equations (37,37) are used in Eq. (29).

Our numerical results are shown in Fig. 5. For simplicity, we have used here the constant $`\mathrm{\Gamma }_\mathrm{\Delta }=120\mathrm{M}\mathrm{e}\mathrm{V}.`$ We notice considerable dependence on $`𝐪`$, as well as different behavior for the transverse and longitudinal cases. The transverse width decreases with $`|𝐪|`$, while the longitudinal does not. At lower values of $`M`$ and $`|𝐪|`$ around $`0.5`$GeV the longitudinal width develops a hill, absent in the transverse case.

The quantity which enters the formula for the dilepton production (see Sect. 5) is the spectral function of the transverse and longitudinal $`\rho `$ mesons, defined as555The presence of $`\sqrt{M^2+𝐪^2}`$ here is related to the presence of $`1/q_0`$ in Eq.(29), i.e. to the fact that we are using widths viewed from the rest frame of the medium. Had we used widths viewed from the $`\rho `$ rest frame, we would have $`1/M`$ in Eq.(29), and factors of $`M`$ instead of $`\sqrt{M^2+𝐪^2}`$ in Eq. (38). Of course, in both cases the resulting spectral functions are equal.

$$𝒜_P=\frac{1}{\pi }\frac{\sqrt{M^2+𝐪^2}\mathrm{\Gamma }_P}{(M^2m_\rho ^2)^2+(M^2+𝐪^2)\mathrm{\Gamma }_P^2},P=T,L,$$
(38)
where $`m_\rho `$ is the position of the pole. The results for $`𝒜_T`$ and $`𝒜_L`$, with $`m_\rho =776`$MeV, are plotted in Fig. 6. We note that while at $`𝐪=0`$ we obviously have $`𝒜_T=𝒜_L`$, at larger values of $`𝐪`$ and at $`M`$ around $`m_\rho `$ the transverse spectral strength becomes dominant. At first glance this may seem surprising, since in Fig. 5 we have seen that at higher $`𝐪`$ we have much larger $`\mathrm{\Gamma }_L`$ than $`\mathrm{\Gamma }_T`$. However, the optimum value of $`\mathrm{\Gamma }_P`$, at which $`𝒜_P`$ has a maximum, is $`\mathrm{\Gamma }_P=(M^2m_\rho ^2)/\sqrt{M^2+𝐪^2}`$. Lower, as well as higher values of $`\mathrm{\Gamma }_P`$ lead to a decrease in $`𝒜_P`$. This explains the behavior of Fig. 6. We note that the transverse spectral strength, $`𝒜_T`$, is concentrated along a ridge extending far into the large-$`𝐪`$ region. Thus, a proper description of propagation at finite and large values of $`𝐪`$ is needed for the description of $`\rho `$ mesons in medium. We note that our results are in qualitative agreement with Ref.  (Fig. 14), with the somewhat different behavior of $`𝒜_T`$, which reaches larger values at higher values of $`|𝐪|`$ in our approach.666Note a factor of $`1/\pi `$ difference in our definition of the spectral functions compared to those of Ref. . We also find qualitative similarity to the results of the altogether different model of Ref.  (Figs. 6,7). There is a difference at larger values of $`|𝐪|`$ for $`𝒜_T`$, manifest in the presence of the rim in our Fig. 6. The results for $`𝒜_L`$ are very similar to .

We should stress here that our construction of the $`\rho `$ spectral function accounts only for the process $`\rho \pi \pi `$ and ignores all other possible contributions, such as e.g. from the $`s`$-channel resonances, studied in Refs. , medium modifications of the pions, etc. Such processes should be included in a more complete calculation. Still, the results of Sect. 5 depend only on the shape of the spectral functions of Fig. 6, and not on the physics leading to their form. Since the spectral functions obtained in many other approaches are quite similar to ours, the results obtained below can be viewed as representative to approaches containing the broadening of the $`\rho `$.

## 5 Dilepton production rate

Measurements of the low-mass dilepton spectra  have shown significant excess above yields from the final-state hadron decays. In this context, the properties of vector mesons (especially of the $`\rho `$) in a hadronic environment become of particular interest, since the Vector Meson Dominance Model is commonly used to make the estimates of the dilepton yields from vector-meson decays. There are numerous analysis of the effect in the literature, both in hydrodynamic approaches , and in transport theories . The dilepton-rate formula  from $`\rho `$ meson decays can be written in a manifestly Lorentz covariant way as follows:

$$\frac{dN}{d^4xdM^2}=\frac{d^3q}{(2\pi )^3}\frac{M}{E_q}\mathrm{\Gamma }_{\rho e^+e^{}}𝒜(M,qu,\rho _B(x))f_\rho \left(\frac{qu}{T(x)}\right),$$
(39)
where $`M`$ is the invariant mass of the lepton pair, equal to the mass of the virtual $`\rho `$ meson, $`𝒜=2𝒜_T+𝒜_L`$ is the spectral function including all polarizations, $`x`$ is a space-time point, $`T(x)`$ is the local value of temperature, and $`E_q=\sqrt{M^2+𝐪^2}`$. The quantity $`\mathrm{\Gamma }_{\rho e^+e^{}}`$ in Eq. (39) is the width for the process $`\rho e^+e^{}`$,

$$\mathrm{\Gamma }_{\rho e^+e^{}}=\frac{4\pi \alpha _{QED}^2m_\rho ^4}{3g_\rho ^2M^5}\left(1\frac{4m_e^2}{M^2}\right)^{\frac{1}{2}}\left(M^2+2m_e^2\right),$$
(40)
where $`\alpha _{QED}`$ is the fine structure constant, and $`m_e`$ is the mass of the electron. Finally, the function $`f_\rho `$ in Eq. (39) is the thermal Bose-Einstein distribution of the $`\rho `$ mesons,

$$f_\rho =\left[\text{exp}\left(\frac{qu2\mu _\pi }{T}\right)1\right]^1,$$
(41)
whith $`\mu _\pi `$ denoting the pion chemical potential , incorporated in several works. This quantity in some sense mimics possible deviations of the system from the chemical equilibrium.

In order to describe the problem as realistically as possible, we will include the effects of the expansion of the medium formed in a relativistic heavy-ion collision. The lepton pairs are formed in a fire cylinder which moves as a whole in the lab system with the rapidity $`\alpha _{FC}`$. For symmetric and central collisions $`\alpha _{FC}`$ is a half of the projectile rapidity in the lab. In its own center-of-mass system (CM), the fire cylinder undergoes a hydrodynamic expansion. In the analysis of such a situation, it is convenient to rewrite Eq. (39) in the variables suited to both the kinematics of the emission process and the geometry of the experimental setup. We introduce $`M_{}=\sqrt{M^2+q_{}^2}`$, the transverse mass of the dilepton pair, $`y^{lab}`$, the rapidity of the pair measured in the lab system, $`𝐮_{}`$, the transverse four-velocity of the fluid element producing dileptons, and $`\alpha ^{lab}`$, the rapidity of this fluid element in the lab. With these variables we have

$$qu=M_{}\sqrt{1+u_{}^2}\text{cosh }(y^{lab}\alpha ^{lab})𝐪_{}𝐮_{}.$$
(42)
The velocity of the fluid element in the lab is a relativistic superposition of the velocity of the fire cylinder in the lab and the hydrodynamic flow considered in the CM system. Thus we have

$$\alpha ^{lab}=\alpha +\alpha _{FC}=\text{arctanh }v_{||}+\alpha _{FC},u_{}=\frac{v_{}\text{cosh}(\alpha )}{\sqrt{1v_{}^2\text{cosh}^2(\alpha )}}.$$
(43)
We note that the velocities $`v_{||}`$ and $`v_{}`$ are defined now in the CM system. They depend on time and space coordinates.

Next, we analyze the kinematic constraint of the CERES experiment. The experimental acceptance cuts can be included with help of the function

$`\mathrm{\Phi }(M,y^{lab},q_{})={\displaystyle \frac{d^2p_1d^2p_2𝑑y_1𝑑y_2\varphi \delta (E_qE_{p_1}E_{p_2})\delta ^{(3)}(𝐪𝐩_1𝐩_2)}{d^2p_1d^2p_2𝑑y_1𝑑y_2\delta (E_qE_{p_1}E_{p_2})\delta ^{(3)}(𝐪𝐩_1𝐩_2)}},`$
(44)
where $`𝐩_{1,\mathrm{\hspace{0.17em}2}}`$ are the momenta of the emitted electrons, $`y_{1,\mathrm{\hspace{0.17em}2}}`$ are the electron rapidities, and $`\varphi `$ is a product of step functions which enforces the experimental setup conditions: $`2.1=\eta _{min}<y_{1,\mathrm{\hspace{0.17em}2}}<\eta _{max}=2.65`$, $`p_{1,\mathrm{\hspace{0.17em}2}}^{}>`$ 200 MeV, and $`\theta _{ee}>`$ 35mrad. Due to the smallness of the electron mass, we can assume here that rapidities and pseudorapidities of the electrons are equal. The construction of the function $`\mathrm{\Phi }(M,y^{lab},q_{})`$ requires a numerical calculation of a two-dimensional integral of a function involving a product of step functions, which is very easily accomplished by a Monte Carlo method. With the inclusion of the experimental acceptance cuts, the dilepton production rate is

$$\frac{dN}{d^4xdM\mathrm{\Delta }\eta }=\frac{2M^2}{(\eta _{max}\eta _{min})}\frac{d^2q_{}}{(2\pi )^3}𝑑y^{lab}\mathrm{\Phi }\mathrm{\Gamma }_{\rho e^+e^{}}𝒜f_\rho .$$
(45)
One should stress here the relevance of the inclusion of the kinematic cuts for the obtained results. The function $`\mathrm{\Phi }`$ influences mostly the overall normalization of the cross section, and not so much the dependence on $`M`$.

In order to calculate the dilepton spectrum, one has to assume a model of the hydrodynamic expansion of the fire cylinder. We adopt the fire-cylinder expansion model of Refs. . It is assumed that the system is in thermal equilibrium up to time $`t_{max}`$, when freeze-out occurs, and the velocities depend on space-time in the following way:

$$v_{||}(t,z)=(v_z+a_zt)\frac{z}{z_{max}(t)},v_{}(t,r)=(v_r+a_rt)\frac{r}{r_{max}(t)},$$
(46)
where

$$z_{max}(t)=z_0+v_zt+\frac{1}{2}a_zt^2,r_{max}(t)=r_0+v_rt+\frac{1}{2}a_rt^2$$
(47)
are the boundaries of the system at time $`t`$. The parameters of the expansion are as follows :

$`t_{max}`$ $`=`$ $`11\mathrm{f}\mathrm{m},z_0=4.55\mathrm{fm},r_0=4.6\mathrm{fm},`$ (48)
$`v_z`$ $`=`$ $`0.5,a_z=0.023\mathrm{fm}^1,v_r=0,a_r=0.05\mathrm{fm}^1,`$ (49)
and the time dependences of the temperature and the baryon density are:

$$T(t)=210\mathrm{M}\mathrm{e}\mathrm{V}\mathrm{exp}\left(\frac{t}{18.26\mathrm{fm}}\right),\rho _B(t)=260/V(t),$$
(50)
where $`V(t)=2\pi z_{max}(t)r_{max}^2(t)`$ is the volume of the fire cylinder at time $`t`$. For the time dependence of the pion chemical potential, $`\mu _\pi (t)`$, we assume a linear rise from $`20\mathrm{M}\mathrm{e}\mathrm{V}`$ at $`t=0`$ to $`80\mathrm{M}\mathrm{e}\mathrm{V}`$ at $`t=t_{max}`$ .

Finally, the yield of leptons produced during the expansion is

$$\frac{dN_1}{dM\mathrm{\Delta }\eta }=\underset{0}{\overset{t_{max}}{}}𝑑t\underset{0}{\overset{r_{max}(t)}{}}2\pi r𝑑r\underset{z_{max}(t)}{\overset{z_{max}(t)}{}}𝑑z\left(\frac{dN}{d^4xdM\mathrm{\Delta }\eta }\right),$$
(51)
where $`dN/(d^4xdM\mathrm{\Delta }\eta )`$ is given by Eq. (45) with all the required substitutions.

In addition to the yields of Eq. (51) one usually adds the contribution from vector mesons which remain after freeze-out. This contribution is equal to

$$\frac{dN_2}{dM\mathrm{\Delta }\eta }=\frac{1}{\mathrm{\Gamma }(M)}\underset{0}{\overset{r_{max}(t_{max})}{}}2\pi r𝑑r\underset{z_{max}(t_{max})}{\overset{z_{max}(t_{max})}{}}𝑑z\left(\frac{dN}{d^4xdM\mathrm{\Delta }\eta }\right),$$
(52)
where $`\mathrm{\Gamma }(M)`$ is the full width of the $`\rho `$ meson with virtual mass $`M`$, given by the formula

$$\mathrm{\Gamma }(M)=\frac{g_{\rho \pi \pi }^2}{48\pi M^2}(M^24m_\pi ^2)^{3/2}$$
(53)
with $`g_{\rho \pi \pi }=5.98`$ giving $`\mathrm{\Gamma }(m_\rho )=150`$MeV. The physical interpretation of formula (52) is that all $`\rho `$ mesons that remain after freeze-out, decay with the yield proportional to the number of mesons and the branching ratio to the dilepton channel. Finally, the full contribution is

$$\frac{dN}{dM\mathrm{\Delta }\eta }=\frac{dN_1}{dM\mathrm{\Delta }\eta }+\frac{dN_2}{dM\mathrm{\Delta }\eta }.$$
(54)
Our numerical results are shown in Fig. 7. The solid line shows the yield from the $`\rho `$ decays, Eq. (54), with the vacuum spectral strength of the $`\rho `$, while the dashed curve shows the calculation with our medium-modified spectral strength of Eq. (38). We can see that the medium effects redistribute the dilepton yields from higher to lower values of $`M`$. This is a typical effect of broadening of the $`\rho ,`$ found in previous investigations. Thus the tendency needed to explain the experimental data is correct. We can see by comparing the dashed line to the data that the calculated yield falls an order of magnitude below the data. Note that other processes, not included in our calculation, contribute in the explored region of $`M`$, for instance the Dalitz decays of mesons, or the $`\omega `$ decays. Also, note that we are not including the “cocktail” contribution of decays of hadrons in our comparison, which would not be consistent. We have found that 15-25% of the model yields comes from decays after freeze-out. The effects of the expansion of the fire cylinder enhance the yields by a few percent. We stress that the overall normalization of the calculated curve is sensitive to the time-integrated volume of the fire cylinder, and to the hydrodynamic expansion parameters .

## 6 Summary and discussion

Here are the main conclusions of our investigation: the medium effects on the $`\rho \pi \pi `$ coupling are large, and dominantly come from the process where the $`\mathrm{\Delta }`$ is excited in the intermediate state. At the nuclear saturation density and at physical mass of the $`\rho `$ the value of the coupling is roughly doubled compared to the vacuum value. The increased coupling leads directly to large widths of the $`\rho `$ meson in medium. We have analyzed the resulting spectral functions for the transverse and longitudinal polarizations, with the result that at higher values of the three-momentum with respect to the medium, the transverse spectral function is much larger from the longitudinal one. Finally, we have applied our model to evaluate the dilepton production from the $`\rho `$ decays in relativistic heavy-ion collisions. We confirm the well-known finding that a larger $`\rho `$ width helps to understand the experimental data by redistributing strength from higher to lower invariant dilepton masses. However, the overall norm of the dilepton yields from the $`\rho `$ decays is almost an order of magnitude too small to explain the experimental data. We remark that the results for the dilepton yields from $`\rho `$ decays are only sensitive to the actual shape of the spectral functions and to the hydrodynamic expansion. In that sense they are not sensitive to detailed modelling of the dynamics, as long as the resulting functions $`𝒜_T`$ and $`𝒜_L`$ are similar. Thus our results confirm the statement of Ref. , namely that hydrodynamic models have problems in explaining the dilepton data unless the hydrodynamic evolution is exceedingly long.

For simplicity, we have used on-shell hadronic couplings throughout our studies. Since the nucleon and the $`\mathrm{\Delta }`$ in the hadronic loop of Fig. 1 can be off-shell, additional coupling structures may be present. In addition, sideways form factors could be included for the particles off-shell. Presently, this has not been done, again for simplicity and from the lack of knowledge as to how to introduce and choose these form factors. One could also include the medium modifications of the meson-baryon couplings, a feature advocated, e.g., in the model of Ref. .

The width of the $`\rho `$ meson picks up contributions not only from the pion-loop diagram, included in our work, but also from other processes. In particular, one can include the $`s`$-channel resonances, as studied e.g. in Refs. . Such processes can and should be included in a more complete and realistic calculation. Also the diagrams of Fig. 1 could in principle be supplied with higher resonances, at the expense of having more not well known parameters.

In addition to the effects of the Fermi sea, studied in this paper, vacuum polarization effects may influence the $`\rho \pi \pi `$ coupling. To have an estimate of these effects, we have done a Walecka-type calculation where in the diagrams of Fig. 1 we have included the free nucleon propagators only. We have found a $``$10% increase of the coupling when the nucleon mass is reduced from 939MeV to 700MeV. We have applied Pauli-Villars regulator with the cut-off parameter of 1GeV in order to truncate high momenta in the loop. Similar order of the effect is found when the Nambu–Jona-Lasinio model is used along the lines of Ref. , and the mass of the quark is scaled down as expected from the medium effects. Thus, vacuum polarization effects on the $`\rho \pi \pi `$ coupling are estimated to be less significant than the Fermi-sea effects presented in this paper.

Concerning the methods applied, we have shown in detail how to obtain very simply the leading-density approximation for diagrams involving loops with density-dependent nucleon propagators. We have also demonstrated how to carry a Lorentz-covariant calculation of the transverse and longitudinal spectral functions with help of the formulas of Sect. 4. These methods and formulas, although very straightforward, are not, to our knowledge, commonly know. They can be useful in studies similar to ours.

One of us (WF) thanks Dariusz Miśkowiec for a discussion of the CERES experimental cuts.

## A Rarita-Schwinger spinors

The Rarita-Schwinger spinors are defined as

$$u^\mu (p,s_\mathrm{\Delta })=\underset{\lambda ,s_N}{}1\frac{1}{2}\lambda s|1\frac{1}{2}\frac{3}{2}s_\mathrm{\Delta }e^\mu (p,\lambda )u_\mathrm{\Delta }(p,s),$$
(A.1)
with

$`e^0(p,\lambda )`$ $`=`$ $`{\displaystyle \frac{\stackrel{}{\epsilon }_\lambda \stackrel{}{p}}{m_\mathrm{\Delta }}},e^i(p,\lambda )=\epsilon _\lambda ^i+{\displaystyle \frac{(\stackrel{}{\epsilon }_\lambda \stackrel{}{p})p^i}{m_\mathrm{\Delta }(E_\mathrm{\Delta }+m_\mathrm{\Delta })}},`$ (A.2)
$`E_\mathrm{\Delta }`$ $`=`$ $`\sqrt{\stackrel{}{p}^2+m_\mathrm{\Delta }^2}.`$ (A.3)
The polarization vectors are defined as

$$\stackrel{}{\epsilon }_0=\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right),\stackrel{}{\epsilon }_\pm =\frac{1}{\sqrt{2}}\left(\begin{array}{c}1\\ i\\ 0\end{array}\right),$$
(A.4)
and

$$u_\mathrm{\Delta }(p,s)=\sqrt{\frac{E_\mathrm{\Delta }+m_\mathrm{\Delta }}{2m_\mathrm{\Delta }}}\left(\begin{array}{c}1\\ \frac{\stackrel{}{\sigma }\stackrel{}{p}}{E_\mathrm{\Delta }+m_\mathrm{\Delta }}\end{array}\right)\chi (s),$$
(A.5)
where $`\chi (s)`$ is the two-component spinor. The spinor $`u^\mu (p,s_\mathrm{\Delta })`$ satisfies the conditions $`\gamma _\mu u^\mu (p,s_\mathrm{\Delta })=0`$ and $`p_\mu u^\mu (p,s_\mathrm{\Delta })=0`$.

## B Isospin algebra

The isospin $`\frac{1}{2}\frac{3}{2}`$ transition matrices are defined through the Clebsch-Gordan coefficients as follows: $`\frac{3}{2},I_3|T^\mu |\frac{1}{2},i_3=\frac{1}{2}1i_3\mu |1\frac{1}{2}\frac{3}{2}I_3`$, with $`i_3`$ and $`I_3`$ denoting the isospin of the nucleon and $`\mathrm{\Delta }`$, respectively. In Cartesian basis the explicit form reads

$$T^1=\left(\begin{array}{cc}\frac{1}{\sqrt{2}}& 0\\ 0& \frac{1}{\sqrt{6}}\\ \frac{1}{\sqrt{6}}& 0\\ 0& \frac{1}{\sqrt{2}}\end{array}\right),T^2=i\left(\begin{array}{cc}\frac{1}{\sqrt{2}}& 0\\ 0& \frac{1}{\sqrt{6}}\\ \frac{1}{\sqrt{6}}& 0\\ 0& \frac{1}{\sqrt{2}}\end{array}\right),T^3=\left(\begin{array}{cc}0& 0\\ \sqrt{\frac{2}{3}}& 0\\ 0& \sqrt{\frac{2}{3}}\\ 0& 0\end{array}\right).$$
(B.1)
where the columns are labeled by $`i_3=\frac{1}{2},\frac{1}{2}`$, left to right, and the rows by $`I_3=\frac{3}{2},\frac{1}{2},\frac{1}{2},\frac{3}{2}`$, top to bottom. The following useful relation holds:

$$T^aT^b=\frac{2}{3}\delta ^{ab}\frac{1}{3}\epsilon ^{abc}\tau _c.$$
(B.2)
The couplings of the $`\mathrm{\Delta }`$ to the isovector proceed via the matrix $`T_\mathrm{\Delta }^\mu `$ defined as $`\frac{3}{2},I_3^{}|T_\mathrm{\Delta }^\mu |\frac{3}{2},I_3=\frac{\sqrt{15}}{2}\frac{3}{2}1I_3^{}\mu |\frac{3}{2}\frac{1}{2}\frac{3}{2}I_3`$. Explicitly, we find

$`T_\mathrm{\Delta }^1`$ $`=`$ $`\left(\begin{array}{cccc}0& \frac{\sqrt{3}}{2}& 0& 0\\ \frac{\sqrt{3}}{2}& 0& 1& 0\\ 0& 1& 0& \frac{\sqrt{3}}{2}\\ 0& 0& \frac{\sqrt{3}}{2}& 0\end{array}\right),T_\mathrm{\Delta }^2=i\left(\begin{array}{cccc}0& \frac{\sqrt{3}}{2}& 0& 0\\ \frac{\sqrt{3}}{2}& 0& 1& 0\\ 0& 1& 0& \frac{\sqrt{3}}{2}\\ 0& 0& \frac{\sqrt{3}}{2}& 0\end{array}\right),`$ (B.11)
$`T_\mathrm{\Delta }^3`$ $`=`$ $`\left(\begin{array}{cccc}\frac{3}{2}& 0& 0& 0\\ 0& \frac{1}{2}& 0& 0\\ 0& 0& \frac{1}{2}& 0\\ 0& 0& 0& \frac{3}{2}\end{array}\right),`$ (B.16)
where the columns are labeled by $`I_3^{}=\frac{3}{2},\frac{1}{2},\frac{1}{2},\frac{3}{2}`$, left to right, and the rows by $`I_3=\frac{3}{2},\frac{1}{2},\frac{1}{2},\frac{3}{2}`$, top to bottom. The conventional factor of $`\sqrt{15}/2`$ in the definition ensures that $`T_\mathrm{\Delta }^3`$ simply measures the third component of the isospin of the $`\mathrm{\Delta }.`$ The spin coupling matrices used in non-relativistic calculations, $`S^i`$ and $`S_\mathrm{\Delta }^i`$, are defined analogously and have exactly the same values as $`T^a`$ and $`T_\mathrm{\Delta }^a`$.

Isospin trace factors for diagrams of Fig. 1 can be now readily obtained. They are equal to $`2`$, $`\frac{2}{3}`$, $`\frac{2}{3}`$, $`\frac{5}{3}`$, $`\frac{5}{3}`$, $`2`$, and $`\frac{4}{3}`$ times $`iϵ^{acb}`$, for diagrams (a),(b),…,(g), respectively. Isospin indices $`a`$, $`b`$, and $`c`$ have the assignment as specified on the text above Eq. (1).

## C Amplitudes for $`𝐪=0`$

For the case $`𝐪=0`$ the amplitude corresponding to the diagram $`(i)`$ of Fig. 1 can be written as

$$A_{(i)}^\mu =\frac{1}{8}\rho _BG_{(i)}\frac{N_{(i)}}{D_{(i)}}(2p^\mu q^\mu ),$$
(C.1)
where $`G_{(i)}`$ are products of coupling constants:

$`G_{(a)}`$ $`=`$ $`{\displaystyle \frac{g_A^2g_\rho }{8F_\pi ^2}},G_{(b)}={\displaystyle \frac{\sqrt{2}g_Af^{}f_{\pi N\mathrm{\Delta }}}{2F_\pi m_\pi ^2}},G_{(c)}={\displaystyle \frac{g_\rho f_{\pi N\mathrm{\Delta }}^2}{2m_\pi ^2}},`$
$`G_{(d)}`$ $`=`$ $`{\displaystyle \frac{g_\rho f_{\pi N\mathrm{\Delta }}^{\mathrm{\hspace{0.17em}\hspace{0.17em}2}}}{m_\pi ^2}},G_{(e)}={\displaystyle \frac{3\sqrt{2}f^{}f_{\pi N\mathrm{\Delta }}f_\mathrm{\Delta }}{2m_\pi ^3}},G_{(f)}={\displaystyle \frac{g_\rho g_A^2}{4F_\pi ^2}},`$
$`G_{(g)}`$ $`=`$ $`{\displaystyle \frac{g_\rho f_{\pi N\mathrm{\Delta }}^2}{m_\pi ^2}}.`$ (C.2)
The formulas for $`N_{(i)}`$ and $`D_{(i)}`$ are very long in the general case, which reflects the presence of many terms in the Rarita-Schwinger propagator. Numerators $`N_{(i)}`$ become manageable in the formal case $`\mathrm{\Gamma }_\mathrm{\Delta }=0`$, and $`m_\pi =0`$, where we find

$`N_{(a)}`$ $`=`$ $`8\left(4m_N^4M^2+\kappa _\rho m_N^2M^4\right),`$
$`N_{(b)}`$ $`=`$ $`16m_N^2\left(m_N+m_\mathrm{\Delta }M\right)M^2\left(m_N+m_\mathrm{\Delta }+M\right)`$
$`\left(2m_N^4+m_N^3m_\mathrm{\Delta }m_Nm_{\mathrm{\Delta }}^{}{}_{}{}^{3}+2m_{\mathrm{\Delta }}^{}{}_{}{}^{4}2m_N^2\left(2m_{\mathrm{\Delta }}^{}{}_{}{}^{2}+M^2\right)\right),`$
$`N_{(c)}`$ $`=`$ $`16m_N\left(m_Nm_\mathrm{\Delta }\right)\left(m_N+m_\mathrm{\Delta }\right)`$ (C.5)
$`\left(m_N^3\left(1+\kappa _\rho \right)m_N^2m_\mathrm{\Delta }2m_Nm_{\mathrm{\Delta }}^{}{}_{}{}^{2}+\left(1+\kappa _\rho \right)m_{\mathrm{\Delta }}^{}{}_{}{}^{3}\right)M^2`$
$`4\left(\begin{array}{c}\left(4+\kappa _\rho \right)m_N^4+4\left(1+\kappa _\rho \right)m_N^3m_\mathrm{\Delta }\\ \left(2+5\kappa _\rho \right)m_N^2m_{\mathrm{\Delta }}^{}{}_{}{}^{2}+2\kappa _\rho m_{\mathrm{\Delta }}^{}{}_{}{}^{4}\end{array}\right)M^4+4\kappa _\rho m_N^2M^6,`$
$`N_{(d)}`$ $`=`$ $`20M^2\left(4m_N^3m_\mathrm{\Delta }m_Nm_{\mathrm{\Delta }}^{}{}_{}{}^{3}4m_{\mathrm{\Delta }}^{}{}_{}{}^{4}+m_N^2\left(7m_{\mathrm{\Delta }}^{}{}_{}{}^{2}+M^2\right)\right),`$
$`N_{(e)}`$ $`=`$ $`40M^2\left(\begin{array}{c}\left(m_Nm_\mathrm{\Delta }\right)\left(m_N+m_\mathrm{\Delta }\right)^3\left(\begin{array}{c}2m_N^4+4m_N^3m_\mathrm{\Delta }\\ 2m_Nm_{\mathrm{\Delta }}^{}{}_{}{}^{3}7m_{\mathrm{\Delta }}^{}{}_{}{}^{4}\end{array}\right)\\ \left(m_N+m_\mathrm{\Delta }\right)\left(\begin{array}{c}3m_N^5+7m_N^4m_\mathrm{\Delta }+4m_N^3m_{\mathrm{\Delta }}^{}{}_{}{}^{2}+\\ m_N^2m_{\mathrm{\Delta }}^{}{}_{}{}^{3}2m_Nm_{\mathrm{\Delta }}^{}{}_{}{}^{4}+2m_{\mathrm{\Delta }}^{}{}_{}{}^{5}\end{array}\right)M^2+\\ m_\mathrm{\Delta }\left(m_N+m_\mathrm{\Delta }\right)\left(2m_N^2+m_Nm_\mathrm{\Delta }2m_{\mathrm{\Delta }}^{}{}_{}{}^{2}\right)M^4+m_N^2M^6\end{array}\right),`$ (C.13)
$`N_{(f)}`$ $`=`$ $`0,`$
$`N_{(g)}`$ $`=`$ $`16\left(\begin{array}{c}4\left(m_Nm_\mathrm{\Delta }\right)m_{\mathrm{\Delta }}^{}{}_{}{}^{2}\left(m_N+m_\mathrm{\Delta }\right)^2+\\ m_N\left(m_N^2+2m_Nm_\mathrm{\Delta }m_{\mathrm{\Delta }}^{}{}_{}{}^{2}\right)M^2\end{array}\right).`$ (C.16)
The expressions for $`D_{(i)}`$ for the case $`\mathrm{\Gamma }_\mathrm{\Delta }=0`$ are

$`D_{(a)}`$ $`=`$ $`\left(4m_N^3m_NM^2\right)\left(m_{\pi }^{}{}_{}{}^{4}+m_N^2M^2\right),`$
$`D_{(b)}`$ $`=`$ $`9m_{\mathrm{\Delta }}^{}{}_{}{}^{2}\left(m_N+m_\mathrm{\Delta }M\right)\left(m_N+m_\mathrm{\Delta }+M\right)\left(m_N^2m_{\mathrm{\Delta }}^{}{}_{}{}^{2}+m_{\pi }^{}{}_{}{}^{2}m_NM\right)`$
$`\times \left(m_{\pi }^{}{}_{}{}^{2}+m_NM\right)\left(m_{\pi }^{}{}_{}{}^{2}+m_NM\right)\left(m_{\mathrm{\Delta }}^{}{}_{}{}^{2}+m_{\pi }^{}{}_{}{}^{2}+m_N\left(m_N+M\right)\right),`$
$`D_{(c)}`$ $`=`$ $`9m_{\mathrm{\Delta }}^{}{}_{}{}^{2}\left(4m_N^3+m_NM^2\right)\left(\left(m_N^2m_{\mathrm{\Delta }}^{}{}_{}{}^{2}+m_{\pi }^{}{}_{}{}^{2}\right)^2m_N^2M^2\right),`$
$`D_{(d)}`$ $`=`$ $`27m_{\mathrm{\Delta }}^{}{}_{}{}^{3}\left(\left(m_N^2m_{\mathrm{\Delta }}^{}{}_{}{}^{2}+m_{\pi }^{}{}_{}{}^{2}\right)^2m_N^2M^2\right),`$
$`D_{(e)}`$ $`=`$ $`27m_{\mathrm{\Delta }}^{}{}_{}{}^{4}(m_N+m_\mathrm{\Delta }M)(m_N+m_\mathrm{\Delta }+M)\times `$
$`\left(\left(m_N^2m_{\mathrm{\Delta }}^{}{}_{}{}^{2}+m_{\pi }^{}{}_{}{}^{2}\right)^2m_N^2M^2\right),`$
$`D_{(f)}`$ $`=`$ $`m_{\pi }^{}{}_{}{}^{4}m_N^2M^2,`$
$`D_{(g)}`$ $`=`$ $`9m_{\mathrm{\Delta }}^{}{}_{}{}^{2}\left(\left(m_N^2m_{\mathrm{\Delta }}^{}{}_{}{}^{2}+m_{\pi }^{}{}_{}{}^{2}\right)^2m_N^2M^2\right).`$ (C.17)


text:
# Greenberger-Horne-Zeilinger nonlocality for continuous variable systems

\[

## Abstract

As a development of our previous work, this paper is concerned with the Greenberger-Horne-Zeilinger (GHZ) nonlocality for continuous variable cases. The discussion is based on the introduction of a pseudospin operator, which has the same algebra as the Pauli operator, for each of the $`N`$ modes of a light field. Then the Bell-CHSH (Clauser, Horne, Shimony and Holt) inequality is presented for the $`N`$ modes, each of which has a continuous degree of freedom. Following Mermin’s argument, it is demonstrated that for $`N`$-mode parity-entangled GHZ states (in an infinite-dimensional Hilbert space) of the light field, the contradictions between quantum mechanics and local realism grow exponentially with $`N`$, similarly to the usual $`N`$-spin cases.

PACS number(s): 03.65.Ud, 03.65.Ta

\]

Under the presumption of local realism, Einstein, Podolsky and Rosen (EPR) argued that it seems possible to assign values with certainty to canonically conjugate variables (e.g., position and momentum) of one particle by performing experiments on another particle of an entangled biparticle system . EPR’s observation is definitely in conflict with quantum mechanics. The entangled states of the biparticle system are now known as the EPR states of continuous variables. In 1965 Bell, using Bohm’s EPR states of discrete variables (e.g., spin) , derived the famous Bell’s inequalities , enabling quantitative tests of quantum mechanics versus local realism. So far, many experiments based on Bohm’s EPR states completely confirmed quantum mechanics . Bell’s theorem without inequalities has also been demonstrated for nonmaximally entangled biparticle states  and multiparticle Greenberger-Horne-Zeilinger (GHZ) states . In this case, the contradictions between quantum mechanics and local realistic theories are of a non-statistical nature; the nonlocality can be manifest in a single measurement. Using three-photon polarization-entangled states, Pan et al. carried out an experimental test of GHZ theorem , confirming again the correctness of quantum mechanics. Besides its theoretical interest, research work on quantum nonlocality also gains great impetus from the rapidly developing field of quantum information theory .

Quantum nonlocality for continuous variables has also attracted much attention recently . Several proposals have been put forward to prepare the original EPR states ; EPR paradox in EPR’s original sense has been realized experimentally . The experiment exploited the nondegenerate optical parametric amplifier (NOPA), which generates the two-mode squeezed vacuum states (“regularized” EPR states)

$$|\mathrm{NOPA}=\frac{1}{\mathrm{cosh}r}\underset{n=0}{\overset{\mathrm{}}{}}(\mathrm{tanh}r)^n|n,n,$$
(1)
with $`r>0`$ denoting the squeezing parameter and $`|nn|n_1|n_2`$ the twin-photon states. In the limit of infinite squeezing, the states $`|\mathrm{NOPA}`$ reduce to a normalized version of the original EPR states . Since the Wigner function  of $`|\mathrm{NOPA}`$ is positive everywhere, it has been argued that the original EPR states may allow a hidden variable description and thus will not exhibit nonlocality . This point of view was dramatically changed by Banaszek and Wódkiewicz (BW) . Within a phase space formulation, BW used the “Bell operator”  based on the joint parity measurements and demonstrated the violations of the Bell inequality due to Clauser, Horne, Shimony and Holt (CHSH) . GHZ nonlocality for continuous variables has also been analyzed , relying on the BW formalism. On the experimental aspect, a recent experiment  observed, for the first time, the violations of the Bell-CHSH inequality by the regularized EPR states produced in a pulsed NOPA within the homodyne detection scheme.

In a recent paper , we have generalized Bell’s inequalities to the continuous variable cases for biparticle systems. Our formalism is based on a new Bell operator, which is a direct analogy of its discrete variable counterpart . In such a physically appealing manner, Bell’s inequalities for both discrete variable (in a two-dimensional Hilbert space) and continuous variable (in an infinite-dimensional Hilbert space) cases take a mathematically similar form. In this paper, we extend the work to entangled multiparty systems. In particular, we will present GHZ theorem for the parity-entangled GHZ states of a three-mode light field. Moreover, the Bell-CHSH inequality for multiparty systems of discrete variables  will be generalized to continuous variable cases.

Let us start with considering a light field with three independent modes $`1`$, $`2`$ and $`3`$. For each mode of the light field, we can introduce the “parity spin” operator (i.e., pseudospin operator acting upon the parity space of photons)

$`\widehat{𝐬}_j`$ $`=`$ $`(s_{jx},s_{jy},s_{jz}),(j=1,2,3)`$ (2)
$`s_{jx}`$ $`=`$ $`s_{j+}+s_j,s_y=i(s_{j+}s_j).`$ (3)
Here

$`s_{jz}`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left[|2n_{jj}2n||2n+1_{jj}2n+1|\right],`$ (4)
$`s_{j+}`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}|2n_{jj}2n+1|=(s_j)^{}`$ (5)
with $`|n_j`$ denoting the usual Fock states of the $`j`$-mode light field. After a simple algebra, it can be shown that

$`[s_{iz},s_{j\pm }]=\pm 2\delta _{ij}s_{i\pm },[s_{i+},s_j]=\delta _{ij}s_{iz}.`$ (6)
$`s_{i\mu }s_{j\nu }+s_{i\mu }s_{j\nu }=2\delta _{ij}\delta _{\mu \nu }.(\mu ,\nu =x,y,z)`$ (7)
The commutation relations in Eq. (7) are identical to their counterparts of the spin-$`1/2`$ systems. Thus each parity spin operator $`\widehat{𝐬}_j`$ acts as the counterpart of the Pauli operator $`\sigma `$. In particular the following relation is valid

$`U_j(\zeta _j,\widehat{𝐧}_j)`$ $`=`$ $`\mathrm{exp}\left(i{\displaystyle \frac{\zeta _j}{2}}\widehat{𝐧}_j\widehat{𝐬}_j\right)`$ (8)
$`=`$ $`\mathrm{cos}{\displaystyle \frac{\zeta _j}{2}}i\widehat{𝐧}_j\widehat{𝐬}_j\mathrm{sin}{\displaystyle \frac{\zeta _j}{2}}`$ (9)
for $`\widehat{𝐬}_j`$ as for $`\sigma `$. Here $`\widehat{𝐧}_j`$ is an arbitrary unit vector. $`U_j(\zeta _j,\widehat{𝐧}_j)`$ are, in fact, the rotation operators in the parity space of photons.

Now one can introduce states of definite parity for each mode

$$|+_j\underset{n=0}{\overset{\mathrm{}}{}}𝒜_n^{(j)}|2n_j,|_j\underset{n=0}{\overset{\mathrm{}}{}}𝒜_n^{(j)}|2n+1_j,$$
(10)
where $`_{n=0}^{\mathrm{}}\left|𝒜_n^{(j)}\right|^2=1`$. Some important properties of these parity states are

$`s_{jz}|\pm _j`$ $`=`$ $`\pm |\pm _j,`$ (11)
$`s_{j\pm }|\pm _j`$ $`=`$ $`0,s_{j\pm }|_j=|\pm _j.`$ (12)
Hence $`|\pm _j`$ are the eigenstates of the parity operators $`s_{jz}`$, with the corresponding eigenvalues $`\pm 1`$; the operators $`s_{j\pm }`$ “flip” the parity of $`|_j`$ and annihilate $`|\pm _j`$. In this aspect, $`|\pm _j`$ are similar to the spin-up and spin-down states of the spin-$`1/2`$ systems. Using the rotation operators as in Eq. (9), one can construct eigenstates of $`\widehat{𝐧}_j\widehat{𝐬}_j`$ for each mode. Now the parity-entangled GHZ state of the light field can be directly constructed as

$$|\mathrm{GHZ}_3=\frac{1}{\sqrt{2}}\left(|+_1|+_2|+_3|_1|_2|_3\right).$$
(13)
Note that this state is defined in an infinite-dimensional Hilbert space.

Following GHZ’s argument , it is directly verifiable that the GHZ state (13) satisfies the eigenvalue equations as

$`s_{1x}s_{2x}s_{3x}|\mathrm{GHZ}_3`$ $`=`$ $`|\mathrm{GHZ}_3,`$ (14)
$`s_{1x}s_{2y}s_{3y}|\mathrm{GHZ}_3`$ $`=`$ $`|\mathrm{GHZ}_3,`$ (15)
$`s_{1y}s_{2x}s_{3y}|\mathrm{GHZ}_3`$ $`=`$ $`|\mathrm{GHZ}_3,`$ (16)
$`s_{1y}s_{2y}s_{3x}|\mathrm{GHZ}_3`$ $`=`$ $`|\mathrm{GHZ}_3.`$ (17)
Thus quantum mechanically, one can determine any operator ($`s_{1x}`$, $`s_{1y}`$, $`s_{2x}`$,$`\mathrm{}`$) by performing appropriate measurements on the other two modes of the light field. To establish a local realistic interpretation of the quantum-mechanical results (17), one assumes that the individual value of any operator ($`s_{1x}`$, $`s_{1y}`$, $`s_{2x}`$,$`\mathrm{}`$) is predetermined, e.g., by local hidden variables, regardless of any set of parity spin measurements on the three modes of the light field. These predetermined values are denoted by $`p_{1x}`$, $`p_{1y}`$, $`p_{2x}`$,$`\mathrm{}`$, with $`p_{jx},p_{jy}=\pm 1`$. To be consistent with Eq. (17), local realistic theories impose the following relations

$`p_{1x}p_{2x}p_{3x}`$ $`=`$ $`1,p_{1x}p_{2y}p_{3y}=+1,`$ (18)
$`p_{1y}p_{2x}p_{3y}`$ $`=`$ $`+1,p_{1y}p_{2y}p_{3x}=+1.`$ (19)
But in fact, these relations are not mutually consistent. For example, the latter three relations in Eq. (19) give $`p_{1x}p_{2x}p_{3x}p_{1y}^2p_{2y}^2p_{2y}^2=p_{1x}p_{2x}p_{3x}=+1`$, which conflicts with the first relation in Eq. (19). Here we have used the fact that $`p_{jy}^2=1`$. Thus quantum-mechanical predictions are incompatible with local realistic theories for the parity-entangled GHZ states (13). Also, the contradictions between quantum mechanics and local realism occur for perfect correlations (i.e., for non-statistical predictions) of the three-mode light field. It is worth pointing out that the GHZ paradox associated with position-momentum variables has also been considered in terms of the Weyl algebra in Ref. .

In a real experiment, the perfect correlation condition in the GHZ argument is very difficult to implement practically. To face this difficulty, one has to rely on the multiparty generalization of Bell’s inequalities, which has been derived originally for discrete variables . Here we generalize the multiparty Bell-CHSH inequality to the present continuous variable case.

First, let us define the Bell operator in the case of a two-mode light field :

$`_2`$ $`=`$ $`(𝐚_1\widehat{𝐬}_1)[(𝐚_2\widehat{𝐬}_2)+(𝐚_2^{}\widehat{𝐬}_2)]`$ (21)
$`+(𝐚_1^{}\widehat{𝐬}_1)[(𝐚_2\widehat{𝐬}_2)(𝐚_2^{}\widehat{𝐬}_2)],`$
where $`𝐚_{1,2}`$ and $`𝐚_{1,2}^{}`$ are four unit vectors. The commutation relations in Eq. (7) lead to, e.g.,

$$(𝐚_1\widehat{𝐬}_1)^2=I_1,$$
(22)
where $`I_1`$ is the unit operator for the mode-$`1`$ of the light field. Equation (22) implies that the eigenvalues of the Hermitian operator $`𝐚_1\widehat{𝐬}_1`$ are $`\pm 1`$. Now the $`N`$-mode Bell operator can be recursively defined by

$`_N`$ $`=`$ $`_{N1}{\displaystyle \frac{1}{2}}[(𝐚_N\widehat{𝐬}_N)+(𝐚_N^{}\widehat{𝐬}_N)]`$ (24)
$`+_{N1}^{}{\displaystyle \frac{1}{2}}[(𝐚_N\widehat{𝐬}_N)(𝐚_N^{}\widehat{𝐬}_N)],`$
where $`\widehat{𝐬}_m`$ ($`m=1`$, $`2`$, $`\mathrm{}`$, $`N`$) are the parity spin operators as defined in Eq. (3), $`𝐚_m`$ and $`𝐚_m^{}`$ are all unit vectors. $`_N^{}`$ defined in Eq. (24) has the same expression as $`_N`$, but with $`𝐚_m𝐚_m^{}`$. Similarly to the discrete variable case , local realistic theories must satisfy the following $`N`$-party Bell-CHSH inequality

$$\left|_N\right|2,$$
(25)
where $`_N`$ is the expectation value of $`_N`$ with respect to any $`N`$-mode state of the light field.

Meanwhile, an upper bound of $`\left|_N\right|`$ can be determined quantum mechanically. By analogy to the $`N`$-party Bell-CHSH inequality for discrete variables , we can prove, using the properties of $`\widehat{𝐬}_m`$ \[see Eq. (7)\], that

$`_N^2`$ $`=`$ $`_{N1}^2{\displaystyle \frac{1}{2}}(1+𝐚_N𝐚_N^{})`$ (28)
$`+_{N1}^2{\displaystyle \frac{1}{2}}(1𝐚_N𝐚_N^{})`$
$`+[_{N1},_{N1}^{}]{\displaystyle \frac{i}{2}}(𝐚_N\times 𝐚_N^{})\widehat{𝐬}_N`$
$``$ $`(1+\left|𝐚_N\times 𝐚_N^{}\right|)\left|_{N1}^2\right|2^{N+1}.`$ (29)
In the last step of Eq. (29), we have used the fact that $`_2^22^{2+1}`$ , from which the upper bound $`2^{N+1}`$ of $`_N^2`$ can be recursively obtained. Therefore quantum mechanics gives an upper bound $`2^{(N+1)/2}`$ of $`\left|_N\right|`$; $`N`$-mode entangled states of the light field can violate the $`N`$-mode Bell-CHSH inequality (25) by a maximal factor of $`2^{(N1)/2}`$, similarly to the discrete variable case . For the two-mode Bell-CHSH inequality, the upper bound $`2^{3/2}`$ of $`\left|_N\right|`$ is known as the Cirel’son bound . Interestingly, the $`N`$-party Bell-CHSH inequality (25), though expressed in terms of discrete variable operators (i.e., the parity spin operators), can be exploited to uncover the GHZ nonlocality of continuous variable states, e.g., the multi-mode entangled states considered in Ref. .

Can an $`N`$-mode entangled state of the light field maximally violate the $`N`$-mode Bell-CHSH inequality (25)? To answer this question, we follow Mermin’s reasoning  and take, as an example, the following $`N`$-mode parity-entangled GHZ state (instead of the multi-mode entangled states as in Ref. )

$$|\mathrm{GHZ}_N=\frac{1}{\sqrt{2}}\left(|+_1|+_2\mathrm{}|+_N|_1|_2\mathrm{}|_N\right).$$
(30)
From Eq. (12) we have

$`(s_{mx}\pm is_{my})|_m`$ $`=`$ $`2|\pm _m,`$ (31)
$`(s_{mx}\pm is_{my})|\pm _m`$ $`=`$ $`0.`$ (32)
Then it can be easily verified that $`|\mathrm{GHZ}_N`$ is an eigenstate of the following operator

$`A`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[{\displaystyle \underset{m=1}{\overset{N}{}}}(s_{mx}+is_{my})+{\displaystyle \underset{m=1}{\overset{N}{}}}(s_{mx}is_{my})\right],`$ (34)
$`A|\mathrm{GHZ}_N=2^{N1}|\mathrm{GHZ}_N`$
with eigenvalue $`2^{N1}`$. Hence, quantum mechanics predicts that

$$\left|{}_{N}{}^{}\mathrm{GHZ}\left|A\right|\mathrm{GHZ}_{N}^{}\right|=2^{N1}.$$
(35)
By contrast, local realistic theories predict a much lower bound on the measured results of $`A`$ for large $`N`$. Under the notion of local realism, there are predetermined values $`p_{1x}`$, $`p_{1y}`$, $`p_{2x}`$,$`\mathrm{}`$ ($`p_{jx},p_{jy}=\pm 1)`$ corresponding to operators $`s_{1x}`$, $`s_{1y}`$, $`s_{2x}`$,$`\mathrm{}`$. Each of these predetermined values of these operators does not depend on other measurements. Consequently, local realistic prediction on the upper bound of $`A`$ reads

$$\left|\mathrm{Re}\underset{m=1}{\overset{N}{}}(p_{mx}+ip_{my})\right|\left|\underset{m=1}{\overset{N}{}}(\pm 1\pm i)\right|=2^{N/2}.$$
(36)
Compared with the quantum-mechanical prediction (35), local realism predicts a bound of $`A`$ lower by a large factor $`2^{N/21}`$, growing exponentially with $`N`$.

As comparison, within the BW formalism the contradictions between quantum mechanics and local realism also grow with increasing number of parties for multi-mode entangled states generated from squeezed states and beam splitters . But the growth seems to decrease for larger number of parties, instead of the exponential growth. In this aspect, the merit of the present formalism is manifest.

To summarize, we have discussed the GHZ nonlocality for continuous variables, as a natural development of our previous work . The introduction of the parity spin operators satisfying the same algebra as the Pauli operator is crucial to our discussion. Then the $`N`$-party Bell operator of continuous variables was constructed in terms of the parity spin operators of the $`N`$-mode light field, resulting in a generalization of the multiparty Bell-CHSH inequality to the continuous variable cases. Following Mermin’s argument, it has been shown that for $`N`$-mode parity-entangled GHZ states (in an infinite-dimensional Hilbert space) of the light field, the contradictions between quantum mechanics and local realism grow exponentially with $`N`$, similarly to the usual $`N`$-spin cases. Within our formulation, there is a striking similarity between the multiparty Bell theorem for both discrete and continuous variables. Our result shows that only partial (i.e., parity spin) information of the whole Hilbert space of continuous variables is sufficient for the mere purpose of uncovering quantum nonlocality. On an experimental aspect, however, the problem of realizing the parity-entangled GHZ states remains open.

Note added.—After the completion of the work, we became aware of a related paper by Massar and Pironio . These authors considered the GHZ paradox for continuous variables in terms of modular and binary variables.

We thank Rob Clifton for bringing Ref.  into our attention. This work was supported by the National Natural Science Foundation of China under Grants No. 10104014, No. 19975043 and No. 10028406 and by the Chinese Academy of Sciences.


text:
# Binary fluid demixing: The crossover region

## I Introduction

Spinodal decomposition occurs when a binary mixture, which forms a homogeneous phase at high temperature, undergoes spontaneous demixing after imposing a sudden quench to below the spinodal temperature of the mixture. Following the quench, after a short period of interdiffusion, the mixture will form domains of different composition, separated by sharply defined interfaces. In its late stage, the local compositions of the fluid domains correspond to those of the two bulk phases at coexistence, while the interfacial tension approaches its equilibrium value. The shape and evolution of the domains will depend on the initial composition of the mixture, as well as on the different physical parameters characterizing the fluids. For the case of a symmetric mixture, in which the amounts and properties of the two immiscible species are the same, then starting from a randomly mixed initial state, bicontinuous structures develop.

As emphasized by Siggiasiggia  and Furukawafurukawa , the physics of spinodal decomposition involves capillary forces, viscous forces and fluid inertia. Local interfacial curvature generates stress which drives fluid motion; this will propagate by viscous force or inertial motion (or both) depending on the parameters of the fluid. If we assume no other physics determines the spinodal process, the parameters determining the fluid behavior are then the interfacial tension, $`\sigma `$, fluid mass density, $`\rho `$, and viscosity $`\eta `$ (for a deep quench, diffusion does not contribute significantly beyond a short initial transient). Any degree of asymmetry in the mixture’s composition or in the dynamical properties of the fluids will provide additional control parameters.

For the completely symmetric case, from the three relevant fluid parameters only one length, $`L_0`$, and one time, $`T_0`$, can be constructed: $`L_0=\eta ^2/(\rho \sigma )`$ and $`T_0=\eta ^3/(\rho \sigma ^3)`$. The absence of other physical mechanisms controlling the spinodal decomposition in the late-stage domain growth leads to the dynamical scaling hypothesissiggia ; furukawa . According to it, if we express the domain size, $`L(T)`$ as a function of time $`T`$ in reduced units using the characteristic length and time scales, the corresponding dimensionless length $`lL/L_0`$, expressed in terms of the dimensionless time $`tT/T_0`$, will be a universal scaling function

$$l=l(t)$$
(1)
the same for all fully symmetric, binary incompressible fluid mixtures. The use of scaling concepts has been of extreme help in understanding and rationalizing the physics underlying complex phenomena. In the present problem, the dynamical scaling states that all fluids will have a unique domain growth rate governed only by the value of $`L/L_0`$, which controls the dominance of viscous or inertial fluid transport. Therefore, this scaling expression provides a systematic and simple way to analyze the spinodal decomposition of any fluid, in terms of its basic physical mechanisms.

Furukawa (following Siggia) furukawa ; siggia  analyzed the limiting behavior of the proposed universal function $`l(t)`$. For small enough $`t`$ fluid inertia is negligible compared to viscous forces, while for large $`t`$ the opposite is trueinertial . Balancing the dominant terms in the Navier-Stokes equation in each case leads to the asymptotes

$$l(t)\{\begin{array}{cc}bt\hfill & ;tt^{}\hfill \\ ct^{2/3}\hfill & ;tt^{}\hfill \end{array}$$
(2)
where $`t^{}`$ stands for the crossover time, when the viscous and inertial contributions are comparable. This can be defined more precisely as the time at which the extrapolated viscous and inertial asymptotes of the $`l(t)`$ cross each other on a log-log plot. For intermediate times around $`t^{}`$, the universal function $`l(t)`$ will exhibit a smooth crossover from the viscous to the inertial asymptotic behavior. It is important to remark that dynamical scaling implies the universality, not just of the asymptotic power laws for the scaling function $`l(t)`$, but of the entire curve, and accordingly of the amplitudes $`b`$, $`c`$ and $`t^{}`$.

Although there is good experimental evidence for the viscous growth regimeexpt1 , the inertial regime is more difficult to study experimentallyviv2prl ; vivjfm . Therefore, simulation techniques have so far been the basic tools to analyse the different growth regimes and validate the scaling hypothesis vivjfm ; simulations ; vivprl ; jury . (The ease with which fully symmetric immiscible fluids can be created on a computer, compared to the laboratory, is also a factor.) Such scaling has allowed classification and comparison of (previously scattered) simulation results as well as a critical analysis of the derived parameters (such as $`b`$), revealing a number of discrepancies in the published data (mainly related to residual diffusion and/or finite size effects)vivjfm . In our view, the existence of a viscous scaling regime (with a universal $`b`$) is now fairly clearly established in three dimensions. (The case of two dimensions is quite different wagner ; wagner2D .) The search for an asymptotic inertial regime has likewise been extended carefully up to Reynolds numbers of order 300 (with Reynolds number defined as $`Re=L\dot{L}\rho /\eta =\dot{l}l`$). Although there is still a controversy regarding the true final asymptote at sufficiently high Reynolds numbergrant ; viv2prl , previous workvivjfm ; vivprl  does establish the existence of a well defined asymptotic inertial regime, at least as far as the $`l(t)`$ curve is concerned, at relatively modest $`Re100`$. The attainment of asymptotic behaviour in this region of $`l(t)`$ (roughly, $`3\times 10^3<l<3\times 10^5`$) is less obvious when one examines velocity statisticsvivjfm  or interfacial dynamicsnjp . However, further exploration at still higher $`Re`$ (say, $`l10^6`$) remains out of reach with present simulation methods.

Compared to the status of the viscous and inertial regimes, the existence of dynamical scaling in the crossover regime remains to be established. So far it was reported in vivjfm ; vivprl  that the crossover region spans four decades in dimensionless time, and that the crossover time $`t^{}`$, although formally of order one, is about $`10^4`$. But the lattice Boltzmann data presented in vivjfm  is in fact relatively sparse for the crossover region, and does not allow a direct check of the proposed collapse onto a universal curve.

There are, however, some results found using dissipative particle dynamics jury  which lie within the crossover regime as determined in vivjfm . Jury et al jury  reported incomplete collapse of this data: each dataset had a similar slope but the curves did not join up on the $`l(t)`$ plot. Their data was found to be broadly compatible with logarithmic deviations from dynamic scaling. Although the more recent (lattice Boltzmann) results obtained in vivprl  suggest that those deviations may instead be due to finite size effects, the validity of dynamical scaling within the crossover region is far from completely established.

Jury et al. noted that the dynamical scaling hypothesis relies on the assumption that viscous forces and inertia are the only physical mechanisms controlling domain growth. Other physics is, however, involved too: we have already mentioned diffusion, which can however be kept negligible, on the domain scale, for deep quenches (see Ref.vivjfm  for a discussion of this). Another mechanism that must be present during the domain growth arises during topological reconnection or pinch-off events. These events are necessary for the coarsening of the bicontinuous structure. They correspond to the contraction of a fluid neck to zero width in a finite time (see e.g. Ref.njp  for visualizations of the evolution and rupture of necks in the different regimes of spinodal decomposition). Unless diffusion intervenes first, the final stage of such a reconnection event always involves a microscopic, molecular, length scale. Moreover, recent studies indicate that such reconnection events follow a distinct scaling of the dynamicseggers ; brenner  that might, in principle, interfere somehow with the scaling law for the domain size. It is not clear whether this could happen only in the crossover region and not in the viscous or inertial asymptotes, but nonetheless, there is enough uncertainty to warrant a detailed study of the crossover domain, which we now present.

## II Model

Spinodal decomposition with a deep quench is essentially a problem of deterministic, isothermal fluid motion coupled to moving interfaces, in which it is important to deal with pinch-off events consistently. As already mentioned, diffusion is irrelevant, and thermal fluctuations can also be disregarded at long times. To study this process we have therefore used the lattice Boltzmann method. This is a simulation technique in which the interfaces emerge as a result of the imposed free energy. This fact has some advantages with respect to other standard numerical techniques to study Navier-Stokes equations. The details of the model for the case of a binary mixture we have used are described elsewherejccpc ; vivjfm ; swift . The model simulates a binary mixture with the model free energy

$$F=𝑑𝐫\left\{\frac{A}{2}\varphi ^2+\frac{B}{4}\varphi ^4+\frac{1}{3}\rho \mathrm{log}\rho +\frac{\kappa }{2}|\varphi |^2\right\},$$
(3)
in which $`A`$, $`B`$ are parameters that control the phase behavior of the mixture. When $`A<0`$ and $`B>0`$ two phases coexist in equilibrium. In Eq.3, $`\varphi `$ is the usual order parameter (the normalized difference in number density of the two fluid species),and if $`A=B`$, the order parameter in the coexisting low temperature phases takes values $`\varphi =\pm 1`$. (This is a matter of convention, not of physics.) The parameter $`\kappa `$ is related to the energy cost of generating a spatial gradient. It determines the value of the interfacial width, $`\xi =5\sqrt{\kappa /2A}`$, and the surface tension $`\sigma =\sqrt{8\kappa A^3/9B^2}`$. Finally, $`\rho `$ is the total fluid density, which remains essentially constant during phase separation. (This is ensured by working at very low Mach numbers.)

The lattice Boltzmann model for a binary mixture involves two velocity distribution functions, $`f`$ and $`g`$, on a discrete lattice with a discrete time relaxational dynamics. The properties of the model are fixed through the equilibrium distributions towards which $`f`$ and $`g`$ relax; the zeroth, first and second moments of $`f`$ determine the density, the momentum flux, and the stress tensor, while moments of $`g`$ determine $`\varphi `$ and the order parameter flux. In this way the dynamics of the distribution functions is connected to the macroscopic behavior of the fluid. It has been shown that the Navier-Stokes and advective-diffusion equations are recovered as the hydrodynamic limit of the corresponding lattice Boltzmann dynamicsswift , at least in the incompressible limitvivjfm . The form of the pressure tensor depends on the free energy model chosenswift .

The spontaneous emergence of interfaces as a result of the imposed free energy, eq.(3), is an appealing feature of this method. However, as a result the interfaces are not structureless, but have finite width $`\xi `$. In practice, the parameters of the free energy are chosen such that $`\xi `$ remains around 3 lattice spacings in all runs; this minimizes anisotropy effects due to the underlying lattice. The finite width imposes restrictions on the length scales that can accurately be studied with this method. In real fluids (deeply quenched) the interfacial width is a few atomic diameters, whereas we are often in the regime of $`\xi L_0`$ (especially for the most inertial runs), which is much larger. For spinodal decomposition this should not matter as long as the dynamical scaling hypothesis holds. However, if the interfacial width turns out to be dynamically relevant, as could possibly apply during pinch-off, the variation between runs of $`\xi /L_0`$ has to be taken into account.

Since we are interested in testing the dynamical scaling hypothesis in the absence of diffusion, the mobility must also be chosen to ensure that, after a short initial transient which is always dominated by diffusion (when the initial domains are formed) diffusion gives an irrelevant contribution to the growth rate. This can be done using protocols described invivjfm .

## III Method and results

We have performed a number of new lattice Boltzmann simulations on large ($`256^3`$) lattices, with parameters chosen to explore the crossover region, as well as a couple of runs deeper into the viscous and inertial regimes to clarify the convergence to the corresponding asymptotic regimes. These add significantly more data to those presented previously in vivprl ; vivjfm , where the focus was on the asymptotes rather than the crossover regime. We start with random initial conditions (effectively infinite temperature), and choose the parameters in the free energy, eq.(3), such that $`B>0`$ and $`A=B`$. This is equivalent to performing a sudden deep quench below the spinodal temperature of the mixture (which corresponds to $`A=0`$). We then follow the subsequent evolution of the system by computing the characteristic domain size $`L(T)`$, which is then converted to reduced physical units.

There are several different ways to measure the domain size in a binary mixturevivjfm . Here we choose to employ a characteristic length scale extracted from the curvature of the interface (see Ref.wagner ). The latter is a tensor, defined by

$$D_{\alpha \beta }=\frac{_{lattice}_\alpha \varphi _\beta \varphi }{_{lattice}\varphi ^2}$$
(4)
where $`\alpha `$ and $`\beta `$ refer to the three spatial coordinates, and $`_\alpha `$ means the spatial derivative with respect to the $`\alpha `$ component of the spatial coordinate. The three eigenvalues of this matrix $`\lambda _i`$ are then three characteristic inverse lengths. From them, the domain size, $`L`$, can be estimated as

$$L=\frac{3}{\lambda _1+\lambda _2+\lambda _3}$$
(5)
Strictly speaking, this is a square of a length; it is related to the product of two characteristic lengths. However, the second relevant length for the order parameter gradients in the system is the interfacial width. This is the same for all the simulations; hence, it contributes a constant factor of order unity to eq.(5). We have monitored the values of the domain size as a function of time after the quench, $`L(T)`$. Before seeking data collapse by converting to reduced physical units, it is very important to select the subset of reliable data from each simulation run. Early in each run one always has a period of interdiffusion, followed by a regime where diffusion and hydrodynamics are both present. Only when the diffusive contribution to the coarsening rate becomes small can the $`L(T)`$ data be expected to scale onto the universal $`l(t)`$ curve. (We removed the early time data by eye but the cutoff values we used were comparable to those that V. Kendon obtained by a more rigorous procedure vivjfm .) On the other hand, it is found that when the domain size is much larger than about a quarter of the simulation box, finite size effects start to be relevant. Hence, from each simulation only a central portion of the overall data can be relied upon; this data covers at most one decade or so in time vivjfm . As a result, the full crossover region cannot be spanned by a single run but only by varying simulation parameters to access $`l(t)`$ one section at a time. Since other parameters, such as $`\xi /L_0`$, are varying during this procedure, there is no guarantee, in general, of a smooth joining up of the curves. In table 1 we show the set of fluid parameters we have used to explore the crossover regime. We have also added a couple of sets that lie in the asymptotic regimes to check that the universal viscous and inertial scaling regimes are indeed recovered.

The initial period of diffusive growth, in each run, means that the reduced physical time must be defined as $`t=(TT_{int})/T_0`$, where the offset time $`T_{int}`$ is not known a priori, and can vary from one run to another. Previously, Kendon et al. developed a careful procedure for fitting this parametervivjfm , separately for each run, using the data from that run only. Although this is an objective procedure, the fitted values of $`T_{int}`$ are very sensitive to the details of the fitting, and the precise location of the resulting $`l(t)`$ data is similarly sensitive. In the current work, therefore, we leave $`T_{int}`$ as a free parameter for each run, and then globally optimize the choices of all the $`T_{int}`$’s to achieve the best data collapse on the $`l(t)`$ curve. This is best done by human eye using an interactive graphics routine. The $`T_{int}`$ values are displayed in table 2; given the uncertainties, these are broadly consistent with the ones obtained previously.

We show in Fig.1 the scaling curve we have obtained in this fashion for all the runs of table 1. One can clearly see that all the simulation data can thereby be made to fall on a universal scaling curve, with overlapping datasets extending from the viscous regime deep into the inertial one. We have also drawn the two lines that correspond to the viscous and inertial asymptotic laws to display how these are attained at either end of the crossover. Since we are using a definition for the domain size different from that of Ref.vivjfm , we have to recompute the amplitudes of the asymptotic law. These are $`b=0.065`$ and $`c=1.3`$, which gives a crossover time $`t^{}=10^4`$. Our crossover region now spans about two decades, somewhat narrower than the previous estimates; in particular, runs that were previously taken to lie in the late crossover (vivjfm, ; njp, , Run29,Run30) are found to be much closer to the asymptote than previously suggested. The reduced crossover width is more in line with normal expectations and suggests that uncertainty in the determination of $`T_{int}`$ by the previous method vivprl ; vivjfm  led to an exaggeration of the crossover.

In Fig.2 we zoom in on the initial part of the crossover region. It is apparent how (to within small errors), the different runs lie on top of each other, as pieces of a continuous universal curve. It is interesting to look at Run13. Although in general, one run spans only one dynamical regime (due to the reasons explained at the beginning of the section), Run13 clearly exhibits a linear growth followed by a deviation, at the beginning of the crossover region. That this is not a finite-size artifact is clear, since runs that lie further in the crossover region coincide with the final portion of Run13.

In Fig.3 we have focused on the late part of the crossover region. Again, one can clearly see how different runs lie on top of each other. The relaxation towards the final, inertial growth law is faster than previously predicted. There is no evidence of any kind of deviations. This nice collapse into a single curve rules out the logarithmic deviations jury  that could be attributed to the emergence of a new length scale in this dynamical regime.

Our procedure for choosing $`T_{int}`$ (by optimizing data collapse) is less objective than the one used previously by Kendon et al.vivjfm  and it is not surprising that the resulting scaling is better. However, the sensitivity to $`T_{int}`$, and the quality of the collapse when this is allowed to float, leads us to conclude that there is no firm evidence against dynamical scaling in the crossover region, despite the previous results of Jury et al.jury .

Of course, since the fitting is done by human eye, there is a danger, perhaps, of seeing scaling where none exists. To check this, we have tried the same procedure in two dimension where scaling is indeed held not to exist (at least, not for bicontinuous morphologies)wagner ; wagner2D . The same approach to the data does not give similar collapse of the two dimensional data, which is reassuring.

As a further validation of our approach, we have also checked our results by computing the domain size using an independent domain size measurement. Specifically, for Runs 3, 4, 7, 13, 20 and 29 we have computed the spherically averaged order parameter structure factor

$$S(k,T)_{|𝐤|=k}\varphi (𝐤,T)\varphi (𝐤,T)d^3𝐤$$
(6)
from it, we can obtain a measure of the mean domain size as the inverse of its first moment,

$$L_\varphi (T)=\frac{S(k,T)𝑑k}{kS(k,T)𝑑k}$$
(7)
In this way we again sample the universal curve from the viscous up to the inertial regime. We have used the values of $`T_{int}`$ obtained from $`L`$ to construct the universal curve for $`L_\varphi `$. In fig. 4 we show the different runs after scaling appropriately, which shows that we are able to collapse all the data on a single universal curve. It exhibits the two asymptotic viscous and inertial regimes. The fact that we are able to use the obtained values of $`T_{int}`$ to construct the universal curve using independent measures of the domain lengths gives confidence that the scaling procedure we have followed is robust. If we compare the $`l(t)`$ curves obtained from the two domain lengths, we see that the ones derived from $`L(T)`$ tend to lie above those corresponding to $`L_\varphi `$. However, in fig.5 we compare the scaled curves obtained using the two length measures for a few runs. For run 29 we have estimated the error in the measure as the mean square deviation of three length-scales measured as $`1/\lambda _1`$, $`1/\lambda _2`$ and $`1/\lambda _3`$ (see eq. (5)), which we display as upper and lower bars. The error in the determination of the domain size, which is around $`5\%`$, does not allow to distinguish, in the late portion of the runs, between the two measures. This analysis also shows that the slight deviations from scaling visible in Fig. 2 are well within the error tolerances expectednote .

## IV Conclusions

In this paper we have analyzed in detail the crossover region of the spinodal decomposition of a binary mixture. By performing a number of new large lattice Boltzmann simulations, we have covered the regime where viscous forces and inertial transport are simultaneously relevant. Throughout the crossover, the data is compatible with the existence of a universal scaling curve $`l(t)`$ in reduced physical units. The crossover region is somewhat narrower than previously reported. This can be attributed to the uncertainties in the behavior in the late part of the crossover arising from the ambiguities in the determination of the fitting parameter $`T_{int}`$. The existence of the scaling curve argues against the suggestion jury  that the physics of pinch-off events can strongly interfere with the domain growth in this regime of parameters.

###### Acknowledgements.

This work was funded in part by EPSRC Grant GR/M56234 and EC Access to Research Infrastructure Contract HPRI-1999-CT-00026 TRACS programme at EPCC. .