text:
# 1 Introduction

## 1 Introduction

In the present article we investigate the implications of the recently published data from the salt phase of measurements of the SNO solar neutrino experiment  for the predictions of the effective Majorana mass in neutrinoless double beta decay, $`\left|<m>\right|`$(see, e.g., ). We update our earlier similar analysis in , which was inspired by the 2002 SNO data . As is well known, the neutrinoless double beta ($`(\beta \beta )_{0\nu }`$) decay, $`(A,Z)(A,Z+2)+e^{}+e^{}`$, is allowed if neutrino mixing, involving the electron neutrino $`\nu _e`$, is present in the weak charged lepton current and the neutrinos with definite mass are Majorana particles (see, e.g., ). Strong evidences for neutrino mixing, i.e., for oscillations of solar electron neutrinos $`\nu _e`$ driven by nonzero neutrino masses and neutrino mixing , have been obtained in the solar neutrino experiments : the pioneering Davis et al. (Homestake) experiment  and in Kamiokande, SAGE, GALLEX/GNO and Super-Kamiokande. These evidences have been spectacularly reinforced during the last two years by the data from the SNO solar neutrino and KamLAND reactor antineutrino experiments. Under the assumption of CPT-invariance, the observed disappearance of reactor $`\overline{\nu }_e`$ in the KamLAND experiment  confirmed the interpretation of the solar neutrino data in terms of $`\nu _e\nu _{\mu ,\tau }`$ oscillations, induced by nonzero neutrino masses and nontrivial neutrino mixing. The KamLAND results practically established the large mixing angle (LMA) MSW solution as unique solution of the solar neutrino problem.

The combined 2-neutrino oscillation analysis of the solar neutrino and KamLAND data identified two distinct solution sub-regions within the LMA solution region - LMA-I,II (see, e.g., ). The best fit values of the two-neutrino oscillation parameters - the solar neutrino mixing angle $`\theta _{}`$ and the mass squared difference $`\mathrm{\Delta }m_{}^2`$, in the two sub-regions - LMA-I and LMA-II, read (see, e.g., ): $`\mathrm{\Delta }m_{}^2^I=7.3\times 10^5\mathrm{eV}^2`$, $`\mathrm{\Delta }m_{}^2^{II}=1.5\times 10^4\mathrm{eV}^2`$, and $`\mathrm{tan}^2\theta _{}^I=\mathrm{tan}^2\theta _{}^{II}=0.46`$. The LMA-I solution was preferred statistically by the data. At 90% C.L. it was found in, e.g., :

$$\mathrm{\Delta }m_{}^2(5.617)\times 10^5\mathrm{eV}^2,\mathrm{tan}^2\theta _{}(0.320.72).$$
(1)
Very recently the SNO collaboration published data from the salt phase of the experiment . For the ratio of the CC and NC event rates, in particular, the collaboration finds: $`R_{CC/NC}=0.306\pm 0.026\pm 0.024`$. Adding the statistical and the systematic errors in quadratures one gets at 3 s.d.: $`R_{CC/NC}0.41`$. As was shown in , an upper limit of $`R_{CC/NC}<0.5`$ implies a significant upper limit on $`\mathrm{\Delta }m_{}^2`$ smaller than $`2\times 10^4\mathrm{eV}^2`$. In the case of interest, $`R_{CC/NC}0.41`$, one finds using the results from : $`\mathrm{\Delta }m_{}^2<1.5\times 10^4\mathrm{eV}^2`$. Thus, the new SNO data on $`R_{CC/NC}`$ implies stringent constraints on the LMA-II solution. A combined statistical analysis of the data from the solar neutrino and KamLAND experiments, including the latest SNO results, shows  that the LMA-II solution is allowed only at 99.13% C.L. Furthermore, the data have substantially reduced the maximal allowed value of $`\mathrm{sin}^2\theta _{}`$, excluding the possibility of maximal mixing 222 The future data from SNO on the day-night effect and the spectrum of $`e^{}`$ from the CC reaction, and the future high statistics data from KamLAND, in principle, can resolve completely the LMA-I - LMA-II solution ambiguity and can constrain further the solar neutrino mixing angle (see  and the references quoted therein). at 5.4 s.d. This has very important implications for the predictions of $`\left|<m>\right|`$.

A 3-$`\nu `$ oscillation analyzes of all available data, including the latest SNO results, were performed in  333Combined 2-neutrino oscillation analyzes of the solar neutrino and KamLAND data were completed earlier in refs. , and in .. It was found, in particular, in  that at 90% C.L.,

$`0.23<\mathrm{sin}^2\theta _{}<0.38\text{for}\mathrm{sin}^2\theta =0.0,`$ (2)
$`0.25<\mathrm{sin}^2\theta _{}<0.36\text{for}\mathrm{sin}^2\theta =0.04,`$ (3)
where $`\theta `$ is the angle limited by the CHOOZ and Palo Verde experiments . The best fit (BF) values in both cases read $`\mathrm{sin}^2\theta _{}=0.30`$. The allowed values of $`\mathrm{\Delta }m_{}^2`$in the LMA-I region have not changed considerably and are given at 90% C.L. by :

$`5.6\times 10^5\mathrm{eV}^2<\mathrm{\Delta }m_{}^2<9.2\times 10^5\mathrm{eV}^2\text{for}\mathrm{sin}^2\theta =0.0,`$ (4)
$`6.1\times 10^5\mathrm{eV}^2<\mathrm{\Delta }m_{}^2<8.5\times 10^5\mathrm{eV}^2\text{for}\mathrm{sin}^2\theta =0.04.`$ (5)
The $`\mathrm{\Delta }m_{}^2`$ best fit value is practically the same for the two values of $`\mathrm{sin}^2\theta `$: $`\mathrm{\Delta }m_{}^27.2\times 10^5\mathrm{eV}^2`$.

There are also strong evidences for oscillations of the atmospheric $`\nu _\mu `$ ($`\overline{\nu }_\mu `$) from the observed Zenith angle dependence of the multi-GeV $`\mu `$-like events in the SuperKamiokande experiment . The combined analysis  of atmospheric neutrino data and the data from the K2K long base line accelerator experiment , shows that at 90% C.L. the neutrino mass squared difference responsible for the atmospheric neutrino oscillations lies in the interval

$$2.0\times 10^3\mathrm{eV}^2<\mathrm{\Delta }m_\mathrm{A}^2<3.2\times 10^3\mathrm{eV}^2.$$
(6)
The $`\mathrm{\Delta }m_\mathrm{A}^2`$ best fit value found in  reads: $`\mathrm{\Delta }m_\mathrm{A}^2|_{\mathrm{BF}}=2.6\times 10^3\mathrm{eV}`$. Let us note that the preliminary results of an improved analysis of the SK atmospheric neutrino data, performed recently by the SK collaboration, gave

$$1.3\times \mathrm{\hspace{0.17em}10}^3\text{eV}^2<|\mathrm{\Delta }m_\mathrm{A}^2|<\mathrm{\hspace{0.17em}3.1}\times \mathrm{\hspace{0.17em}10}^3\text{eV}^2,90\%\mathrm{C}.\mathrm{L}.,$$
(7)
with best fit value $`|\mathrm{\Delta }m_\mathrm{A}^2|=2.0\times 10^3`$ eV2. Adding the K2K data , the authors  find the same best fit value and

$$1.55\times \mathrm{\hspace{0.17em}10}^3\text{eV}^2<|\mathrm{\Delta }m_\mathrm{A}^2|<\mathrm{\hspace{0.17em}2.60}\times \mathrm{\hspace{0.17em}10}^3\text{eV}^2,90\%\mathrm{C}.\mathrm{L}.$$
(8)
The last parameter relevant for our further discussion is the neutrino mixing angle $`\theta `$, limited by CHOOZ and Palo Verde experiments. The precise limit on $`\theta `$ is $`\mathrm{\Delta }m_\mathrm{A}^2`$ dependent (see, e.g, ). For the values of $`\mathrm{\Delta }m_\mathrm{A}^2`$found in  (see eq. (6)), the upper bound on $`\mathrm{sin}^2\theta `$ at 90% (99.73%) C.L. reads: $`\mathrm{sin}^2\theta <0.03(0.05)`$. Using the latest SK preliminary result, one gets at 90% (99.73%) C.L. from a combined 3-$`\nu `$ oscillation analysis of the solar neutrino, CHOOZ and KamLAND data :

$$\mathrm{sin}^2\theta <0.047(0.074).$$
(9)
Under the assumptions of 3-neutrino mixing, for which we have compelling evidences from the experiments with solar and atmospheric neutrinos and from the KamLAND experiment, of massive Majorana neutrinos and of $`(\beta \beta )_{0\nu }`$-decay generated only by the (V-A) charged current weak interaction via the exchange of the three Majorana neutrinos $`\nu _j`$, the effective Majorana mass in $`(\beta \beta )_{0\nu }`$-decay of interest is given by (see, e.g., ):

$$\left|<m>\right|=\left|m_1|U_{\mathrm{e1}}|^2+m_2|U_{\mathrm{e2}}|^2e^{i\alpha _{21}}+m_3|U_{\mathrm{e3}}|^2e^{i\alpha _{31}}\right|,$$
(10)
where $`U_{ej}`$, $`j=1,2,3`$, are the elements of the first row of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) neutrino mixing matrix  $`U`$, $`m_j>0`$ is the mass of the Majorana neutrino $`\nu _j`$, and $`\alpha _{21}`$ and $`\alpha _{31}`$ are two Majorana CP-violating phases . One can express  the masses $`m_{2,3}`$ and the elements $`U_{ej}`$ respectively in terms of $`m_1`$, $`\mathrm{\Delta }m_{}^2`$ $`\mathrm{\Delta }m_\mathrm{A}^2`$, and of $`\theta _{}`$, $`\theta `$ (see further).

In ref.  we have analyzed the implications of the results of the solar neutrino experiments, including the 2002 SNO data, which favored the LMA MSW solution of the solar neutrino problem with $`\mathrm{tan}^2\theta _{}<1`$, for the predictions of $`\left|<m>\right|`$. Neutrino mass spectra with normal mass hierarchy, with inverted hierarchy and of quasi-degenerate type were considered (see, e.g., ). From the fact that $`\mathrm{cos}2\theta _{}0.26`$, which followed (at 99.73% C.L.) from the analysis of the solar neutrino data performed in , we found significant lower bounds on $`\left|<m>\right|`$in the cases of neutrino mass spectrum of quasi-degenerate and inverted hierarchical type: $`\left|<m>\right|>0.03`$ eV and $`\left|<m>\right|>8.5\times 10^3`$ eV, respectively. We have also found that if the neutrino mass spectrum were hierarchical (with inverted hierarchy), $`\left|<m>\right|`$were limited from above: $`\left|<m>\right|<8.2\times 10^3(8.0\times 10^2)`$ eV. These results led us to conclude that a measured value of $`\left|<m>\right|>10^2`$ eV could provide fundamental information on the type of the neutrino mass spectrum. It was also shown that such a result could provide a significant upper limit on the mass of the lightest neutrino $`m_1`$ as well. In refs.  the possibilities to determine the type of neutrino mass spectrum and to get information of CP-violation in the lepton sector, associated with Majorana neutrinos from a measurement of $`\left|<m>\right|`$have been studied in greater detail. In these articles, in particular, the uncertainty in the determination of $`\left|<m>\right|`$due to an imprecise knowledge of the relevant nuclear matrix elements and prospective experimental errors in the measured values of the neutrino oscillation parameters, entering into the expression of $`\left|<m>\right|`$, were taken into account. In  the results of analyzes of the solar and atmospheric neutrino and the first KamLAND data were used in obtaining predictions for $`\left|<m>\right|`$444A concise discussion of the relevant formalism, of the status of the searches for $`(\beta \beta )_{0\nu }`$-decay and of the physics potential of the $`(\beta \beta )_{0\nu }`$-decay experiments as well as of the predictions for $`\left|<m>\right|`$before the publication of the latest SNO data, can be found in ..

In the present Addendum, we update the predictions for $`\left|<m>\right|`$derived in ref.  and the conclusion reached in the indicated article by taking into account the implications of the recently announced data from the salt phase measurements of the SNO experiment. We also comment very briefly how the results obtained in  are modified.

## 2 Predictions for the Effective Majorana Mass Parameter $`\left|<m>\right|`$

The predicted value of $`\left|<m>\right|`$ depends in the case of $`3\nu `$ mixing on the oscillation parameters $`\mathrm{\Delta }m_\mathrm{A}^2`$, $`\theta _{}`$, $`\mathrm{\Delta }m_{}^2`$ and $`\theta `$ (see, e.g., ). Following , we will use the convention in which $`m_1<m_2<m_3`$. In this convention one has $`\mathrm{\Delta }m_\mathrm{A}^2\mathrm{\Delta }m_{31}^2`$. In what regards $`\mathrm{\Delta }m_{}^2`$, there are two possibilities: $`\mathrm{\Delta }m_{}^2\mathrm{\Delta }m_{21}^2`$, or $`\mathrm{\Delta }m_{}^2\mathrm{\Delta }m_{32}^2`$. The former corresponds to neutrino mass spectrum with normal hierarchy, while the latter - to neutrino mass spectrum with inverted hierarchy. Obviously, $`m_3=\sqrt{m_1^2+\mathrm{\Delta }m_\mathrm{A}^2}`$. If $`\mathrm{\Delta }m_{}^2\mathrm{\Delta }m_{21}^2`$, one has $`m_2=\sqrt{m_1^2+\mathrm{\Delta }m_{}^2}`$ and the following relations are valid: $`|U_{\mathrm{e1}}|^2=\mathrm{cos}^2\theta _{}(1|U_{\mathrm{e3}}|^2)`$, $`|U_{\mathrm{e2}}|^2=\mathrm{sin}^2\theta _{}(1|U_{\mathrm{e3}}|^2)`$, and $`|U_{\mathrm{e3}}|^2\mathrm{sin}^2\theta `$. In the case of neutrino mass spectrum with inverted hierarchy, we have $`m_2=\sqrt{m_1^2+\mathrm{\Delta }m_\mathrm{A}^2\mathrm{\Delta }m_{}^2}`$ and $`|U_{\mathrm{e2}}|^2=\mathrm{cos}^2\theta _{}(1|U_{\mathrm{e1}}|^2)`$, $`|U_{\mathrm{e3}}|^2=\mathrm{sin}^2\theta _{}(1|U_{\mathrm{e1}}|^2)`$, $`|U_{\mathrm{e1}}|^2\mathrm{sin}^2\theta `$. Given $`\mathrm{\Delta }m_{}^2`$, $`\mathrm{\Delta }m_\mathrm{A}^2`$, $`\theta _{}`$ and $`\mathrm{sin}^2\theta `$, the value of $`\left|<m>\right|`$ depends strongly on the type of the neutrino mass spectrum and on the value of the lightest neutrino mass, $`m_1`$, as well as on the values of the two Majorana CP-violation phases present in the PMNS matrix, $`\alpha _{21}`$ and $`\alpha _{31}`$ (see eq. (10)).

As in ref. , we consider the predictions of $`\left|<m>\right|`$for the following three specific types of neutrino mass spectrum: i) normal hierarchical (NH), corresponding to a spectrum with normal hierarchy and $`m_1m_{2,3}`$, ii) inverted hierarchical (IH), characterized by inverted hierarchy and $`m_1m_{2,3}`$, and iii) quasi-degenerate spectrum (QD) which is realized if $`m_1m_2m_3m_0`$ and $`m_{1,2,3}^2\mathrm{\Delta }m_\mathrm{A}^2,\mathrm{\Delta }m_{}^2`$. Let us note that in the case of the QD spectrum, $`\left|<m>\right|`$ is essentially independent of $`\mathrm{\Delta }m_\mathrm{A}^2`$and $`\mathrm{\Delta }m_{}^2`$, and, as long as the Majorana CP-violation phases $`\alpha _{21}`$ and $`\alpha _{31}`$ are not constrained, the two possibilities, $`\mathrm{\Delta }m_{}^2\mathrm{\Delta }m_{21}^2`$ and $`\mathrm{\Delta }m_{}^2\mathrm{\Delta }m_{32}^2`$, lead effectively to the same predictions for the allowed range of values of $`\left|<m>\right|`$.

In Tables 1 and 2 we give the i) maximal predicted value of $`\left|<m>\right|`$ in the case of NH neutrino mass spectrum, ii) the minimal and maximal values of $`\left|<m>\right|`$ for the IH spectrum, and iii) the minimal value of $`\left|<m>\right|`$ for the QD spectrum. The indicated values of $`\left|<m>\right|`$are calculated for the best fit values and for the 90% C.L. allowed ranges of $`\mathrm{\Delta }m_\mathrm{A}^2`$from  and , and of $`\mathrm{sin}^2\theta _{}`$ and $`\mathrm{\Delta }m_{}^2`$ in the LMA-I solution region, obtained in the analysis in ref. .

In Figs. 1 and 2 we show the allowed ranges of values of $`\left|<m>\right|`$as a function of $`m_1`$ for the cases of spectrum with normal and inverted hierarchy. The predictions for $`\left|<m>\right|`$are obtained by using the best fit (Fig. 1), and the allowed at 90% C.L. (Fig. 2), values of $`\mathrm{\Delta }m_{}^2`$, $`\theta _{}`$ and $`\mathrm{\Delta }m_\mathrm{A}^2`$from refs.  and  and for three fixed values of $`\mathrm{sin}^2\theta `$.

Let us recall that $`\left|<m>\right|_{\mathrm{min}}^{\mathrm{IH},\mathrm{QD}}`$ are given approximately by  (see also ):

$$\left|<m>\right|_{\mathrm{min}}^{\mathrm{IH}}\sqrt{\mathrm{\Delta }m_\mathrm{A}^2}\mathrm{cos}^2\theta |\mathrm{cos}2\theta _{}|,$$
(11)
$$\left|<m>\right|_{\mathrm{min}}^{\mathrm{QD}}m_0\left|\mathrm{cos}^2\theta \mathrm{cos}2\theta _{}\mathrm{sin}^2\theta \right|.$$
(12)
According to , we have at 99.73% C.L.: $`\mathrm{\Delta }m_\mathrm{A}^2>1.1\times 10^3\mathrm{eV}^2`$ and $`\mathrm{sin}^2\theta <0.074`$. The combined analysis of the solar neutrino data, including the latest SNO results, shows  that i) the possibility of $`\mathrm{cos}2\theta _{}=0`$ is excluded at more than 5 s.d., ii) the best fit value of $`\mathrm{cos}2\theta _{}`$ is $`\mathrm{cos}2\theta _{}=0.40`$, and iii) that at 95% C.L. one has for $`\mathrm{sin}^2\theta =0(0.04)`$, $`\mathrm{cos}2\theta _{}>0.22(0.24)`$. These new results firmly establish the existence of significant lower bounds on $`\left|<m>\right|`$in the cases of IH and QD neutrino mass spectra. That has fundamental implications for the searches for $`(\beta \beta )_{0\nu }`$-decay.

A comparison of the results for $`\left|<m>\right|`$, obtained using the best fit values of $`\mathrm{\Delta }m_{}^2`$, $`\theta _{}`$ and $`\mathrm{\Delta }m_\mathrm{A}^2`$, with those reported in refs.  shows that the predictions for $`\left|<m>\right|`$did not change considerably. This is a consequence of the fact that the best fit values of $`\mathrm{\Delta }m_{}^2`$and $`\mathrm{sin}^2\theta _{}`$ are not very different from the values used as input in refs. . With the recent preliminary result $`\mathrm{\Delta }m_\mathrm{A}^2|_{\mathrm{BF}}=2.0\times 10^3\mathrm{eV}^2`$ of the improved analysis of the SK atmospheric neutrino data  one gets somewhat smaller values of $`\left|<m>\right|`$in the case of IH spectrum. For $`\mathrm{sin}^2\theta 0`$, for instance, one finds $`17\times 10^3\mathrm{eV}<\left|<m>\right|^{\mathrm{IH}}<44\times 10^3\mathrm{eV}`$. The ranges of $`\left|<m>\right|`$in the NH and QD spectra depend weakly on $`\mathrm{\Delta }m_\mathrm{A}^2`$in this case 555The effective Majorana mass $`\left|<m>\right|`$for the NH spectrum practically does not depend on $`\mathrm{\Delta }m_\mathrm{A}^2`$if $`\mathrm{sin}^2\theta `$ is sufficiently small, so that $`\sqrt{\mathrm{\Delta }m_{}^2}\mathrm{sin}^2\theta _{}\sqrt{\mathrm{\Delta }m_\mathrm{A}^2}\mathrm{sin}^2\theta `$. Even for $`\mathrm{sin}^2\theta `$ close to its 90% C.L. upper limit of 0.047, the change in $`\left|<m>\right|_{\mathrm{max}}^{\mathrm{NH}}`$ due to the lower value of $`\mathrm{\Delta }m_\mathrm{A}^2|_{\mathrm{BF}}`$ amounts at most to 6 %. and therefore $`\left|<m>\right|_{\mathrm{max}}^{\mathrm{NH}}`$ and $`\left|<m>\right|_{\mathrm{min}}^{\mathrm{QD}}`$ reported in Table 1 are essentially the same as the ones reported in refs. .

According to the 3-$`\nu `$ combined analysis of the solar neutrino and KamLAND and CHOOZ data , for $`\mathrm{sin}^2\theta =0(0.04)`$, the lower bound on $`\mathrm{cos}2\theta _{}`$ at 90% C.L. reads: $`\mathrm{cos}2\theta _{}|_{_{\mathrm{MIN}}}>0.24(0.28)`$. Therefore the main conclusion that in the case of the IH and the QD spectra there exist significant lower bounds on $`\left|<m>\right|`$ , not only still holds, but is considerably strengthened, as we have already emphasized. More specifically, in the case of IH spectrum we get for $`\mathrm{sin}^2\theta =0(0.04)`$ $`\left|<m>\right|_{\mathrm{min}}^{\mathrm{IH}}=0.010(0.011)\mathrm{eV}`$ if we use eq. (6), and $`\left|<m>\right|_{\mathrm{min}}^{\mathrm{IH}}=0.0087(0.0099)\mathrm{eV}`$ utilizing the preliminary result given in eq. (8).

In the case of QD neutrino mass spectrum, a larger lower bound on $`\mathrm{cos}2\theta _{}`$ implies a larger value of $`\left|<m>\right|_{\mathrm{min}}^{\mathrm{QD}}`$, which now reads: $`\left|<m>\right|_{\mathrm{min}}^{\mathrm{QD}}0.0430.048\mathrm{eV}`$. This should be compared with the value found in : $`\left|<m>\right|_{\mathrm{min}}^{\mathrm{QD}}0.03\mathrm{eV}`$.

Let us note that for the presently allowed at 90% C.L. values of $`\mathrm{\Delta }m_{}^2`$, $`\theta _{}`$ and $`\mathrm{\Delta }m_\mathrm{A}^2`$, and for $`\mathrm{sin}^2\theta >0.03`$, there exists a lower bound on $`\left|<m>\right|`$in the case of NH spectrum, provided $`m_1<<10^3`$ eV: one finds $`\left|<m>\right|_{_{\mathrm{max}}}^{\mathrm{NH}}>\mathrm{few}\times 10^4\mathrm{eV}`$. A complete cancellation of the different terms contributing to $`\left|<m>\right|`$is allowed in this case only if

$$\mathrm{sin}^2\theta >\sqrt{\frac{(\mathrm{\Delta }m_{}^2)_{_{\mathrm{MIN}}}}{(\mathrm{\Delta }m_\mathrm{A}^2)_{_{\mathrm{MAX}}}}}(\mathrm{sin}^2\theta _{})_{_{\mathrm{MIN}}}.$$
(13)
Given the 90% C.L. allowed ranges of the oscillation parameters entering into the above inequality, the “cancellation” condition can be satisfied for $`\mathrm{sin}^2\theta >0.038`$. Lower allowed values of $`\mathrm{\Delta }m_\mathrm{A}^2`$or a more stringent upper limit on $`\mathrm{sin}^2\theta `$ would further strengthen this result. It should be emphasized, however, that one can have $`\left|<m>\right|0`$ even for $`\mathrm{sin}^2\theta <0.030`$ if $`m_1`$ is sufficiently large (see Figs. 1 and 2). More generally, since at present the two Majorana CP-violation phases $`\alpha _{21}`$ and $`\alpha _{31}`$ are not constrained and the existing upper bounds on $`m_1`$ are not sufficiently stringent, one can always have $`\left|<m>\right|=0`$ in the case of neutrino mass spectrum with normal hierarchy .

The 95% C.L. allowed ranges of $`\mathrm{\Delta }m_{}^2`$and $`\theta _{}`$ differ only marginally from those derived at 90% C.L., while the upper limit on $`\mathrm{sin}^2\theta `$ changes from 0.047 to 0.053 . The intervals of allowed values of $`\mathrm{\Delta }m_\mathrm{A}^2`$at 95% C.L., according to the older and the latest analyzes  and  do not differ considerably from the 90% C.L. ones, eqs. (6) and (8), and are given respectively by $`\mathrm{\Delta }m_\mathrm{A}^2=(1.83.4)\times 10^3\mathrm{eV}^2`$ and $`\mathrm{\Delta }m_\mathrm{A}^2=(1.42.8)\times 10^3\mathrm{eV}^2`$. These results imply that the predictions for $`\left|<m>\right|`$, obtained using the 95% C.L. allowed ranges of the neutrino oscillation parameters, differ insignificantly from the predictions based on the 90% C.L. ranges of the parameters, which are illustrated in Table 2 and Fig. 2.

Similar conclusion about the existence of a significant and robust lower bound on $`\left|<m>\right|`$, $`\left|<m>\right|>10^2`$ eV, in the case of IH neutrino mass spectrum, has been reached also in , where a $`\chi ^2`$method of analysis was employed.

## 3 Constraining the Type of Neutrino Mass Spectrum and/or CP-violation Associated with Majorana Neutrinos

In refs.  the possibilities to discriminate between the different types of neutrino mass spectrum and to get information about CP-violation associated with Majorana neutrinos by means of a measurement of $`\left|<m>\right|`$were studied. The uncertainty in the relevant nuclear matrix elements and prospective experimental errors in the values of the oscillation parameters, in $`\left|<m>\right|`$, and for the case of QD spectrum - in $`m_1`$, were taken into account. Here we update the results of  using the same notations. We denote by $`(\left|<m>\right|_{\mathrm{exp}})_{_{\mathrm{MIN}}}`$the value of $`\left|<m>\right|`$obtained from a measurement of the $`(\beta \beta )_{0\nu }`$-decay half life-time of a given nucleus by using the largest physical nuclear matrix element, and by $`\zeta `$ the “theoretical uncertainty” in $`\left|<m>\right|`$due to the imprecise knowledge of the nuclear matrix element. Thus, an experiment measuring a $`(\beta \beta )_{0\nu }`$-decay half-life time will determine a range of values of $`\left|<m>\right|`$corresponding to

$$(\left|<m>\right|_{\mathrm{exp}})_{_{\mathrm{MIN}}}\mathrm{\Delta }\left|<m>\right|\zeta ((\left|<m>\right|_{\mathrm{exp}})_{_{\mathrm{MIN}}}+\mathrm{\Delta }),$$
(14)
where $`\mathrm{\Delta }`$ denotes the experimental error in the measurement of $`\left|<m>\right|`$. A part of the analysis in  was performed by using the best fit values of the solar and the atmospheric neutrino oscillation parameters and including 1 s.d. (3 s.d.) uncertainties of 5% (15%) on $`\mathrm{tan}^2\theta _{}`$ and $`\mathrm{\Delta }m_{}^2`$, and of 10% (30%) on $`\mathrm{\Delta }m_\mathrm{A}^2`$. We follow  in the choice for the ranges of the indicated parameters. In Table 3 we report the predicted i) maximal value of $`\left|<m>\right|`$ in the case of NH spectrum, ii) the minimal and maximal values of $`\left|<m>\right|`$ for the IH neutrino mass spectrum, and iii) the minimal value of $`\left|<m>\right|`$ for the QD spectrum.

In order to be possible to discriminate between the NH and the IH spectrum, the following inequality should be fulfilled : $`\zeta \left|<m>\right|_{\mathrm{max}}^{\mathrm{NH}}<\left|<m>\right|_{\mathrm{min}}^{\mathrm{IH}}`$. The latter implies an upper limit on $`\zeta `$. For the currently favored values of the neutrino oscillation parameters and $`\mathrm{sin}^2\theta <0.03`$, the NH spectrum can be distinguished from the IH one even if $`\zeta 3`$. If $`\mathrm{sin}^2\theta >0.03`$, $`\left|<m>\right|_{\mathrm{max}}^{\mathrm{NH}}`$ can be larger, as Table 3 illustrates, and somewhat smaller values of $`\zeta `$ could be required. Similarly, since $`\left|<m>\right|_{\mathrm{min}}^{\mathrm{IH}}\sqrt{\mathrm{\Delta }m_\mathrm{A}^2}|\mathrm{cos}2\theta _{}|`$, a shift of $`\mathrm{\Delta }m_\mathrm{A}^2`$to smaller values would lead to stronger constraints on $`\zeta `$. In the “worst possible case” in which we allow a 3 s.d. error on $`\mathrm{\Delta }m_{}^2`$, $`\mathrm{\Delta }m_\mathrm{A}^2`$and $`\mathrm{tan}^2\theta _{}`$, and use the best fit values of the solar neutrino oscillation parameters, $`\mathrm{\Delta }m_\mathrm{A}^2=2.0\times 10^3\mathrm{eV}^2`$ and $`\mathrm{sin}^2\theta =0.04`$, it is necessary to have $`\zeta <2.3`$.

The possibility to discriminate between the NH and the QD spectra leads to a less stringent condition on $`\zeta `$ than the condition permitting to distinguishing between the NH and the IH spectra: the former is satisfied even for values of $`\zeta `$ exceeding 3. The IH spectrum can be distinguished from the QD spectrum only if $`\zeta <1.5`$ 666 The upper bound on $`\zeta `$ varies from 1.1 to 1.6 according to the different values of $`\mathrm{\Delta }m_\mathrm{A}^2`$and $`\mathrm{sin}^2\theta `$ (see Table 3). The constraint is less stringent for smaller values of $`\mathrm{\Delta }m_\mathrm{A}^2`$and $`\mathrm{sin}^2\theta `$., unless additional information on neutrino masses is provided by the $`^3\mathrm{H}`$ $`\beta `$-decay experiments and/or cosmological and astrophysical measurements.

We update also the conditions on $`\zeta `$ and $`\mathrm{\Delta }`$ which would permit to rule out, or establish, the NH, IH and the QD mass spectrum. The next generation of $`(\beta \beta )_{0\nu }`$-decay experiments are planned to reach a sensitivity of $`\left|<m>\right|0.010.03\mathrm{eV}`$. The QD mass spectrum can be ruled out if $`\zeta <\left|<m>\right|_{\mathrm{min}}^{\mathrm{QD}}/((\left|<m>\right|_{\mathrm{exp}})_{_{\mathrm{MIN}}}+\mathrm{\Delta })`$. For the prospected sensitivity of $`\left|<m>\right|0.01\mathrm{eV}`$ and $`\mathrm{\Delta }0.01\mathrm{eV}`$, the requirement on $`\zeta `$ reads $`\zeta <3`$. In the less favorable case in which $`((\left|<m>\right|_{\mathrm{exp}})_{_{\mathrm{MIN}}}+\mathrm{\Delta })0.05\mathrm{eV}`$, an extremely accurate knowledge of the nuclear matrix elements, i.e. $`\zeta 1`$, would be necessary. In order to establish that the spectrum is of the QD type, the following condition has to be fulfilled: $`((\left|<m>\right|_{\mathrm{exp}})_{_{\mathrm{MIN}}}\mathrm{\Delta })0.2\mathrm{eV}`$, as it implies that $`m_00.2\mathrm{eV}`$. This requirement is satisfied if, e.g., the measured $`(\left|<m>\right|_{\mathrm{exp}})_{_{\mathrm{MIN}}}0.3`$ eV with an experimental error $`\mathrm{\Delta }<0.1`$ eV. The NH spectrum can be excluded provided the measured value of $`\left|<m>\right|`$ is larger than the predicted upper limit on $`\left|<m>\right|`$for this type of spectrum, $`(\left|<m>\right|_{\mathrm{exp}})_{_{\mathrm{MIN}}}\mathrm{\Delta }>\left|<m>\right|_{\mathrm{max}}^{\mathrm{NH}}`$. Given the expected sensitivities on $`\left|<m>\right|`$, this condition can be realized if, e.g., $`(\left|<m>\right|_{\mathrm{exp}})_{_{\mathrm{MIN}}}`$ is at least 3-11 times larger than $`\left|<m>\right|_{\mathrm{max}}^{\mathrm{NH}}`$, i.e. ($`(\left|<m>\right|_{\mathrm{exp}})_{_{\mathrm{MIN}}}0.0150.035`$), and the experimental error amounts to at most $`\mathrm{\Delta }0.010.03\mathrm{eV}`$. A larger measured value of $`(\left|<m>\right|_{\mathrm{exp}})_{_{\mathrm{MIN}}}`$would allow to exclude the NH spectrum even for larger values of $`\mathrm{\Delta }`$. The IH spectrum can be ruled out if $`(\left|<m>\right|_{\mathrm{exp}})_{_{\mathrm{MIN}}}\mathrm{\Delta }>\left|<m>\right|_{\mathrm{max}}^{\mathrm{IH}}`$. For an experimental error on $`\left|<m>\right|`$, $`\mathrm{\Delta }0.01;0.03;0.05\mathrm{eV}`$, this condition is satisfied if $`(\left|<m>\right|_{\mathrm{exp}})_{_{\mathrm{MIN}}}>0.06;0.08;0.1\mathrm{eV}`$. Establishing that the spectrum is of the IH type is quite demanding and requires a measurement of $`\left|<m>\right|`$ with an error $`\mathrm{\Delta }<0.015`$ eV 777Let us note that unless there is additional information on the type of neutrino mass hierarchy or on the range of allowed values of $`m_1`$, it will be in principle impossible to distinguish the case of IH spectrum from the one with partial normal hierarchy  ($`\mathrm{\Delta }m_{}^2\mathrm{\Delta }m_{21}^2`$ and $`0.01\mathrm{eV}<m_1<0.2\mathrm{eV}`$) using only a measurement of $`\left|<m>\right|`$, as the predicted allowed ranges of $`\left|<m>\right|`$in the two cases overlap (see Figs. 1 and 2)..

The possibility of establishing CP- violation in the lepton sector due to Majorana CP-violating phases has been studied in detail in ref. . It was found that it requires quite accurate measurements of $`\left|<m>\right|`$and of $`m_1`$, and holds only for a limited range of values of the relevant parameters. For the IH and the QD spectra, which are of interest for this analysis, the “just CP-violation” region  \- an experimental point in this region would signal unambiguously CP-violation associated with Majorana neutrinos, is larger for smaller values of $`\mathrm{cos}2\theta _{}`$. As the present best fit values of $`\mathrm{\Delta }m_{}^2`$, $`\mathrm{\Delta }m_\mathrm{A}^2`$and especially of $`\mathrm{sin}^2\theta _{}`$ are very similar to those used in , the conclusions reached in ref.  still hold. More specifically, proving that CP-violation associated with Majorana neutrinos takes place requires, in particular, a relative experimental error on the measured value of $`\left|<m>\right|`$not bigger than (15 – 20)%, a “theoretical uncertainty” in the value of $`\left|<m>\right|`$due to an imprecise knowledge of the corresponding nuclear matrix elements smaller than a factor of 2, a value of $`\mathrm{tan}^2\theta _{}>0.55`$, and values of the relevant Majorana CP-violating phases ($`\alpha _{21}`$, $`\alpha _{32}`$) typically within the ranges of $`(\pi /23\pi /4)`$ and $`(5\pi /43\pi /2)`$.

## 4 Conclusions

In the present Addendum, we have updated the predictions for the effective Majorana mass in $`(\beta \beta )_{0\nu }`$-decay $`\left|<m>\right|`$derived in ref.  by taking into account the implications of the recently announced results from the salt phase measurements of the SNO experiment. The combined analyzes of the solar neutrino data, including the latest SNO results, lead to new relatively stringent constraints on the solar neutrino mixing angle $`\theta _{}`$: i) the possibility of $`\mathrm{cos}2\theta _{}=0`$ is excluded at more than 5 s.d., ii) the best fit value of $`\mathrm{cos}2\theta _{}`$ is found to be $`\mathrm{cos}2\theta _{}=0.40`$, and iii) at 95% C.L. one has $`\mathrm{cos}2\theta _{}>0.22`$. These new results firmly establish the existence of significant lower bounds on $`\left|<m>\right|`$in the cases of IH and QD neutrino mass spectra, which in turn has fundamental implications for the searches for $`(\beta \beta )_{0\nu }`$-decay. Using, e.g., the 90% C.L. allowed ranges of the values of the solar and atmospheric neutrino oscillation parameters one finds for the IH and QD spectra, respectively: $`\left|<m>\right|>0.010`$ eV and $`\left|<m>\right|>0.043`$ eV. The lower bounds obtained utilizing the 95% C.L. allowed values of the parameters do not differ substantially from those given above.

We have also updated the earlier results in  on the possibilities i) to discriminate between the different types of neutrino mass spectrum (NH vs IH, NH vs QD and IH vs QD), and ii) to get information about CP-violation induced by the two Majorana CP-violating phases in the PMNS mixing matrix, if a value $`\left|<m>\right|0`$ is measured in the $`(\beta \beta )_{0\nu }`$-decay experiments of the next generation. The remarkable physics potential of the future $`(\beta \beta )_{0\nu }`$-decay experiments for providing quantitative information, in particular, on the type of the neutrino mass spectrum and on the CP-violation associated with Majorana neutrinos, can be fully exploited only if the values of the relevant $`(\beta \beta )_{0\nu }`$decay nuclear matrix elements are known with a sufficiently small uncertainty.

## 5 Acknowledgments

S.T.P. would like to thank D. Vignaud and the members of APC Institute at College de France, Paris, where parts of the work on the present study were done, for kind hospitality. This work was supported in part by the Italian MIUR under the program “Fenomenologia delle Interazioni Fondamentali” (S.T.P.) and by the DOE Grant DE-FG03-91ER40662 and the NASA Grant ATP02-0000-0151 (S.P.).


format: text

text:
# Dynamics in a Bistable-Element-Network with Delayed Coupling and Local Noise

## 1 Introduction

Stochastic rate processes in bi- or multi-stable systems lead to many interesting phenomena observed in various scientific areas ranging from physics to social science, and have thus been studied for a long time.

Here we consider a network of stochastically driven bistable elements whose distinct feature is a time-delayed coupling. The time delays are considered as uniform and the network elements are assumed to be highly interconnected, so that the connectivity can be approximated by a global all to all coupling.

The dynamics of the network is numerically explored by using a Langevin model and analytically by using a Gaussian approximation and a dichotomous model, derived from the corresponding Fokker-Planck equations and Master equation, respectively.

## 2 Langevin model

The Langevin model consists of $`N`$ equations, each describing the overdamped motion of a particle in a bistable potential in the presence of noise and global coupling to a time-delayed mean field $`X(t\tau )=N^1_{i=1}^Nx_i(t\tau )`$,

$$\dot{x}_i(t)=x_i(t)x_i(t)^3+\epsilon X(t\tau )+\sqrt{2D}\xi (t),$$
(1)
where $`\tau `$ is the time delay, $`\epsilon `$ is the coupling strength of the feedback and $`D`$ denotes the variance of the Gaussian fluctuations $`\xi (t)`$.

## 3 Gaussian Approximation

In order to theoretically study the dynamical properties of a globally coupled set of noisy bistable elements (with no time delay), Desai and Zwanzig Desai and Zwanzig (1978) derived a hierarchy of equations for the cumulant moments of the distribution function from the multi-dimensional Fokker-Planck equation for the joint probability distribution for all elements. For large noise intensities, when the statistics of individual elements are approximately Gaussian, this hierarchy can be truncated. Applying this approach to our system yields the following set of equations for the mean field $`X`$ and the variance $`M=N^1(x_iX)^2`$,

$`\dot{X}`$ $`=`$ $`XX^33XM+\epsilon X(t\tau ),`$
$`{\displaystyle \frac{1}{2}}\dot{M}`$ $`=`$ $`M3X^2M3M^2+D.`$
## 4 Dichotomous Model

To study the dynamics of a single bistable element with time-delayed feedback Tsimring and Pikovsky Tsimring and Pikovsky (2001) used a dichotomous approximation. In the dichotomous approximation intra-well fluctuations of $`x_i`$ are neglected, so that a bistable element can be replaced by a discrete two-state system. Applying this complementary approach to our system, allows us to express the mean field dynamics in terms of the hopping rates $`p_{12,21}`$ which denote the probability of a bistable element to change its state form $`1`$ to $`+1`$ and vice versa, respectively. The equation for the mean field reads

$$\dot{X}=p_{12}p_{21}(p_{21}+p_{12})X,$$
(3)
where the hopping probabilities are given by the Kramers transition rate (Kramers, 1940), which in the limit of small noise and small coupling strength reads (cf. Tsimring and Pikovsky (2001))

$$p_{12,21}=\frac{\sqrt{23\epsilon X(t\tau )}}{2\pi }\mathrm{exp}\left(\frac{14\epsilon X(t\tau )}{4D}\right).$$
(4)
## 5 Phase Diagram

A numerical study of the Langevin model (Eq. 1) shows that the system undergoes ordering transitions and demonstrates multistability. That is, for a strong enough positive coupling the system exhibits a non-zero stationary mean field and a variety of stable oscillatory states are accessible for positive and negative feedback. While the transition to the non-zero stationary mean field is second order (continuous), the type of the oscillatory transitions depends on the system parameters and can be first order (discontinuous), associated with hysteretic behavior, or second order.

A linear stability analysis of Eq. (3) yields the critical coupling strengths $`ϵ`$ and the frequencies $`\omega `$ of the accessible oscillatory states. It is found that the oscillation periods can be arbitrary small meaning that a separation of time scales of the macroscopic mean field dynamics and the microscopic hopping dynamics of the constituent bistable elements is not possible.

A comparison of the Langevin model with the Gaussian approximation and the dichotomous theory shows that while near the transition points the system dynamics is strongly non-Gaussian, it is in this regime well described by a two-state model (see Fig. 1), which allows for a complete analysis of the bifurcations of the trivial equilibrium.

## 6 Coherence Resonance

The considered system exhibits the phenomena of coherence resonance and array-enhanced resonance (see Fig. 2).

Both, Kramers random switching frequency $`p`$ (see 4) and the frequency of the oscillatory states $`\omega `$, resulting from the coupling with the time-delayed mean field, depend on the noise strength, i.e., $`p=p(D)`$ and $`\omega =\omega (D)`$. Thus, the noise can be tuned so that the random hoping between the potential wells of the bistable oscillators is synchronized with the periodic modulation of the mean-field. This statistical synchronization takes place when $`\omega =\pi p`$, where the coherence of the oscillatory states becomes maximal.

For the $`N=1`$ case the coherence resonance was observed by Tsimring and Pikovsky (2001). We observe that the resonance phenomenon not only persists in globally coupled networks with large $`N`$, but is enhanced, a property which was found in other systems and is sometimes referred to as array-enhanced resonance (Zhou et al., 2001).

This work was supported by the Swiss National Science Foundation (D.H.) and by the U.S. Department of Energy, Office of Basic Energy Sciences under grant DE-FG-03-96ER14592 (L.T.).


text:
# Numerical simulation of vacuum particle production: applications to cosmology, dynamical Casimir effect and time-dependent non-homogeneous dielectrics

## I Introduction

One of the most intriguing and fascinating features of quantum field theory resides in the non-trivial nature of its vacuum states. Quantum fluctuations present in the vacuum are responsible for non-classical effects that can be experimentally detected. The most well known of such phenomena is perhaps the Casimir effect . In this setting, the vacuum energy between two static metallic plates is shifted due to quantum fluctuations, resulting in an attractive force between the plates. In general non-static situations, changes in a system’s parameters may lead to a redefinition of the natural vacuum state. As a result, a system initially in the vacuum may at later stages display a non-zero particle content. The possibility of this type of particle production mechanism has been identified and studied in a large variety of systems. In a cosmological context, particle creation may occur in back-hole radiation - the Hawking effect , and as a consequence of an expanding gravitational background . In the laboratory it is possible that dynamical generalisations of the Casimir effect may lead to observable effects . Alternatively, there have been efforts to reproduce the Hawking effect using specific matter systems . This is the case of slow-light experiments , black-hole sonic analogues in Bose-Einstein condensates  and dielectric media .

The calculations involved in determining the physical outcome of particle creation processes, though trivial to state, are often hard or impossible to complete. Usually one is required to find the solution of a set of space and/or time-dependent field equations, with initial conditions covering a complete basis of functions. Even when relying on simplifying approximations, the set of problems for which solutions can be found is considerably limited. On the other hand, linear partial differential equations are relatively easy to solve numerically, and the above task is within present computational capabilities. In this paper we explore this possibility and introduce a fully numerical scheme for studying particle production in general settings. We will use examples with known solutions to help developing and testing the required numerical techniques. These will then be applied to new cases, extending original predictions and leading to new results.

In Section II, we rely on a simple cosmological example to introduce notation and discuss concepts relevant for the rest of the paper. We use this system to describe the general numerical approach and we test it, by comparing the results with known analytical predictions. In Section III, we apply these techniques to the case of particle creation by a reflecting moving boundary in a one-dimensional cavity. For an oscillating motion we obtain particle production rates and spectra, for both resonant and non-resonant frequencies, in the short and long-time regimes. In Section IV, we study a system with time dependent optical properties, mimicking a bubble of dielectric material expanding into the vacuum. From the numerical data we determine the dependence of the final particle number on the expansion velocity and strength of the dielectric. Finally, in Section V we discuss our results and suggest further applications. Appendix A contains a description of the numerical algorithm used in Section III to solve the field equations for moving boundary conditions.

## II A First Example: Particle Creation in an Expanding Background

In this section we will look at a simple example of particle creation in a cosmological system. In particular, we will study the evolution of a free complex field $`\varphi `$ in an expanding flat universe. Here and throughout the paper we will follow the conventions and notation of reference.

For a general gravitational background, the equation of motion for a scalar field is given by

$$(g)^{1/2}_\mu [(g)^{1/2}g^{\mu \nu }_\nu \varphi ]+m^2\varphi =0,$$
(1)
where $`g_{\mu \nu }`$ is the space-time metric and $`g=|detg_{\mu \nu }|`$.

The invariant scalar product of any two solutions of Eq. (1) is defined as

$$(\varphi _1,\varphi _2)=i_\mathrm{\Sigma }\left[\varphi _1_\mu \varphi _2^{}_\mu \varphi _1\varphi _2^{}\right]\sqrt{g_\mathrm{\Sigma }}𝑑\mathrm{\Sigma }^\mu ,$$
(2)
where the integral is taken over any space-like surface of the space-time ($`d\mathrm{\Sigma }^\mu `$ is the time-like unit-vector orthogonal to the surface). Using Eq. (2) we can construct orthonormal basis of solutions of the equations of motion, obeying:

$$(u_i,u_j)=\delta _{ij},(u_i^{},u_j^{})=\delta _{ij},(u_i,u_j^{})=0$$
(3)
Expanding the field operator in terms of a given basis we obtain the familiar creation and annihilation operators $`a_i^{}`$, $`a_i`$:

$$\varphi (x)=\underset{i}{}[a_iu_i(x)+a_i^{}u_i^{}(x)]$$
(4)
Particle states are defined in respect to $`a_i^{}`$, $`a_i`$. In general we will be looking at systems which are static and flat for both early and late times - the so called in-out type problems. For these cases we can single-out in a natural way, two different basis of solutions, $`\{u_i^{\mathrm{in}}\}`$ and $`\{u_j^{\mathrm{out}}\}`$ which tend respectively as $`tt_{\mathrm{in}}`$ and $`tt_{\mathrm{out}}`$, to the usual flat-space positive frequency modes. The field can be quantised in terms of any of these two basis, leading to different definitions of creation and annihilation operators and hence to different definitions of particle states. In most cases of interest, the particle content of a given state will change depending on the set of modes we chose as basis. In particular, we will consider the case where the system is in the vacuum state corresponding to the in modes. In respect to the out particle states - the natural choice for late times - this state may have a non-zero particle content. The particle spectrum is given by

$$N_i=\underset{j}{}|\beta _{ij}|^2$$
(5)
where $`N_i`$ is the number of particles in the mode $`i`$. The $`\beta _{ij}`$’s are a subset of the so-called Bogolugov coefficients which relate the two basis and can be evaluated using :

$$\alpha _{ij}=(u_i^{\mathrm{out}},u_j^{\mathrm{in}}),\beta _{ij}=(u_i^{\mathrm{out}},[u_j^{\mathrm{in}}]^{})$$
(6)
In this paper, we will look into how to tackle the calculation of the particle spectrum for general in-out problems, using a fully numerical approach. We will start by evaluating numerically the in modes $`\{u_i^{\mathrm{in}}\}`$. This is done by solving the equation of motion Eq. (1), taking as initial condition the asymptotic, flat-space mode solutions corresponding to the Minkowsky vacuum. These are propagated up to $`t_{\mathrm{out}}`$, when the new stationary flat configuration is reached. The resulting in modes can then be contracted with the out modes, by evaluating numerically the integral Eq. (2) at $`t=t_{\mathrm{out}}`$. In this way we obtain the Bogolugov coefficients Eq. (6). Other physical quantities, such as the particle spectrum can be evaluated in terms of these.

To illustrate the method, we look at a simple example of a $`1+1`$ space-time with a time dependent scale factor. In conformal coordinates, the metric is given by

$$\mathrm{d}s^2=a^2(\eta )(\mathrm{d}\eta ^2\mathrm{d}x^2),$$
(7)
and the field equations of motion are

$$\frac{^2\varphi }{\eta ^2}=\frac{^2\varphi }{x^2}m^2a^2(\eta )\varphi $$
(8)
For this system, the time independent inner product of two solutions is given by

$$(\varphi _1,\varphi _2)=i_{\mathrm{}}^+\mathrm{}\left[\varphi _1_\eta \varphi _2^{}_\eta \varphi _1\varphi _2^{}\right]𝑑x$$
(9)
We take $`a(\eta )a_{\mathrm{in}}`$ as $`\eta \mathrm{}`$, and $`a(\eta )a_{\mathrm{out}}`$ as $`\eta +\mathrm{}`$. In these limits, the corresponding mode basis have the usual flat space-time form:

$`u_k^{\mathrm{in}}{\displaystyle \frac{1}{2\sqrt{\pi \omega _{\mathrm{in}}}}}e^{ikxi\omega _{\mathrm{in}}\eta },\omega _{\mathrm{in}}=\sqrt{k^2+m^2a_{\mathrm{in}}^2}`$ (10)
$`u_k^{\mathrm{out}}{\displaystyle \frac{1}{2\sqrt{\pi \omega _{\mathrm{out}}}}}e^{ikxi\omega _{\mathrm{out}}\eta },\omega _{\mathrm{out}}=\sqrt{k^2+m^2a_{\mathrm{out}}^2}`$ (11)
For the scale-factor $`a(\eta )`$ we choose a particular form for which there is a known analytical solution :

$$a^2(\eta )=\frac{1}{2}(a_{\mathrm{out}}^2+a_{\mathrm{in}}^2)+\frac{1}{2}(a_{\mathrm{out}}^2a_{\mathrm{in}}^2)\mathrm{tanh}(\rho \eta )$$
(12)
Taking Eq. (10) as initial condition, we evolve the equation of motion Eq. (8) up to a time when $`aa_{\mathrm{out}}`$. The numerical solutions are then contracted with the out modes Eq. (11), using Eq. (9) to obtain the Bogolugov coefficients. For our choice of scale-factor these can be evaluated analytically. The derivation can be found in , the result being given by

$`|\alpha _{kk^{}}|^2=\delta _{kk^{}}{\displaystyle \frac{\mathrm{sinh}^2[\pi (\omega _{\mathrm{out}}+\omega _{\mathrm{in}})/(2\rho )]}{\mathrm{sinh}(\pi \omega _{\mathrm{in}}/\rho )\mathrm{sinh}(\pi \omega _{\mathrm{out}}/\rho )}}`$ (13)
$`|\beta _{kk^{}}|^2=\delta _{kk^{}}{\displaystyle \frac{\mathrm{sinh}^2[\pi (\omega _{\mathrm{out}}\omega _{\mathrm{in}})/(2\rho )]}{\mathrm{sinh}(\pi \omega _{\mathrm{in}}/\rho )\mathrm{sinh}(\pi \omega _{\mathrm{out}}/\rho )}}`$ (14)
In Fig. 1 we compare the numerical Bogolugov coefficients with the expected analytical values, for a particular choice of parameters $`m=1`$, $`\rho =10.`$, $`a_{\mathrm{in}}=1.`$ and $`a_{\mathrm{out}}=10.`$ We evolved $`40`$ modes for a period of time $`t_{\mathrm{out}}t_{\mathrm{in}}=0.6`$, enough for the scale factor to change from its initial to final value with good accuracy. The equation of motion Eq. (8) was discretized using a standard Leap-Frog algorithm, in a box of physical size $`L=100`$ with periodic boundary conditions. The lattice spacing and time-step were set to $`dx=0.25`$ and $`dt=0.0001`$. The numerical data matches very well the analytical prediction, confirming the accuracy of the numerical approach for the chosen set of parameters.

Note that only diagonal elements of the Bogolugov matrix are non-zero, a consequence of the invariance of the system under spatial translations. This implies that the particle spectrum coincides with the diagonal of the Bogolugov matrix, i.e. $`N_k=|\beta _{kk}|^2`$. Since Eq. (8) is separable in space and time, we would have the choice of evolving independently each mode’s amplitude in Fourier space. Numerically this approach would have been considerably simpler. Our method is crucial though, in situations where this option is not available. In the next sections we will apply it to a few of such cases, systems that have both explicit space and time-dependency.

## III Particle Production by a Moving Mirror in a Cavity

In this section we will look at particle production by a moving mirror in a closed cavity. This system has attracted the attention of several authors, displaying many complex and non-trivial features despite its apparent simplicity. More importantly, it is the most natural setting for an experimental observation of particle creation from a vacuum state (see  and references therein).

We will consider a scalar field confined to a one-dimensional box with one of the boundaries fixed at $`x=0`$ and the other moving according to a given trajectory $`X(t)`$ (with $`|\dot{X}(t)|<1`$). For the sake of simplicity, and in order to take advantage of conformal invariance we assume the field to be massless, its equation of motion being given by

$$\frac{^2\varphi }{t^2}=\frac{^2\varphi }{x^2}$$
(15)
The boundary conditions for the field are obtained by constraining it to vanish on both walls at all times, i.e. $`\varphi (0,t)=\varphi (X(t),t)=0`$. By doing so we assume the walls to be perfect reflectors.

In the case of a stationary mirror, time translational invariance makes it possible to find a set of mode solutions with positive definite frequency. These are given by

$`u_n(x,t)={\displaystyle \frac{1}{\sqrt{n\pi }}}e^{i\omega _nt}\mathrm{sin}(\omega _nx),`$ (16)
$`\omega _n={\displaystyle \frac{n\pi }{L}}n=1,\mathrm{},\mathrm{}`$ (17)
where $`L`$ is the length of the static box. These modes and their complex conjugates form a complete set of solutions of Eq. (15), orthonormal in respect to the simplectic product

$$(\varphi _1,\varphi _2)=i_0^L\left[\varphi _1_t\varphi _2^{}_t\varphi _1\varphi _2^{}\right]𝑑x$$
(18)
The stationary solution allows us to define a general class of in-out problems. We will consider mirror trajectories such that $`X(t)`$ is a constant $`X_0`$ for $`t<0`$, evolves arbitrarily for a period of time, going back to rest at $`X_f`$ for $`t>t_f`$. The time and space dependence introduced in the system by the motion of the mirror, leads to particle production without the need to consider additional external fields. This phenomenon is known as the dynamical Casimir effect.

In order to determine the final particle content of the box for a given $`X(t)`$, one must find a set of solutions for the equation of motion, using Eq. (27) as initial condition. This problem can be formally solved by taking advantage of the conformal invariance of Eq. (15). G. Moore showed  (see also ) that by changing coordinates, the general problem can be mapped into the stationary case. In this way a complete set of solutions can be found, with the general form $`u_n(x,t)\mathrm{exp}[in\pi R(tx)]\mathrm{exp}[in\pi R(t+x)]`$. $`R(u)`$ is a function to be determined in terms of the particular motion of the mirror, according to

$$R[t+X(t)]=R[tX(t)]+2$$
(19)
This equation is of course non-trivial to solve in generality. One approach is to find $`R(u)`$ in terms of $`X(t)`$ via a perturbative expansion in the velocity of the mirror . This technique is valid only in the adiabatic limit $`\dot{X}<<1`$ and gives no meaningful results in the relativistic regime. Only for a few specific cases of the mirrors’s trajectory have analytical solutions been found .

Using the numerical procedure outlined in Section II the problem can be solved for any given arbitrary mirror trajectory. The known analytical solutions will be useful in providing tests for the algorithms and in gauging the numerical requirements for good accuracy. The general procedure will be as before: the in-modes Eq. (27) are evolved numerically up to the final time $`t_\mathrm{f}`$, when the result is contracted with the basis of stationary modes for the final length $`X_\mathrm{f}`$. The main difficulty lies in obtaining a numerical solution for the wave equation with moving boundaries. Since the size of the box changes with time we are forced to vary the number of points in the simulation lattice. These discontinuous jumps in the lattice size introduce errors in the numerical solution. We overcame the problem by developing an improved Leap-Frog algorithm that smoothes the values of the field derivatives on the boundary - a detailed discussion of the method can be found in Appendix A.

Finally we should stress that although the focus here will be on one dimensional massless field theories, the same methods can easily be extended to higher dimensions and massive fields . In those cases very few results are known  since one cannot rely on conformal invariance to simplify the problem. For those systems, a numerical approach may be the only way to tackle accurately the long time regime.

### A Uniformly Moving Mirror

Here we will look briefly at a simple type of trajectory for which the in-out problem can be solved analytically , namely the case when the mirror moves with constant velocity. The mirror evolves from rest at $`t=0.`$ moving according to $`X(t)=X_0+vt`$ up to final time $`t_f`$, when it becomes stationary again. The function $`R(u)`$ obeying Eq. (19) for this trajectory was fully determined by Castagnino and Ferraro . These authors showed that $`R`$ is a piece-wise linear function with derivative discontinuities. This is a consequence of the discontinuous change in the velocity of the mirror for initial and final times. The final spectrum for created particles was estimated to be of the form $`N(k)v^2/k`$, in the limit of large $`k`$ and velocities much lower than the speed of light.

In Fig. 2 we confirm this estimate. The results were obtained by simulating a box with initial physical length $`X_0=50.`$ expanding with constant speed $`v`$ for a total time $`t_f=50.`$ We evaluated the final particle spectrum for $`20`$ modes, for values of $`v`$ ranging from $`0.1`$ to $`0.9`$. As seen in the figure, the spectrum for high $`k`$ is well fitted by a power law of the form $`N(k)=Bk^\alpha `$, with $`B`$ and $`\alpha `$ depending on $`v`$. As expected, $`\alpha `$ is close to $`1`$, varying between $`1.1`$ for $`v=0.1`$ and $`1.3`$ for $`v=0.9`$. $`B`$ also depends quadratically on the expansion velocity for small $`v`$. When relativistic speeds are reached (about $`v.7`$) this simple relation stops holding, with $`B`$ increasing faster and leading to a stronger rate of particle production.

Note that the $`1/k`$ spectrum leads to a logarithmically divergent total particle number. This divergence is a consequence of the discontinuous changes in the velocity of the mirror . As a rule of thumb we may expect particle production to be related to the mirror’s acceleration. As this diverges for the initial and final times, so does the total particle number. In realistic cases however, the divergence should disappear, as the acceleration remains bounded for all times.

Finally let us note that for a box that expands continuously with constant $`v`$, with no initial or final stationary points, an exact solution for Eq. (19) can easily obtained :

$$R(u)=\frac{1}{\mathrm{arctan}v}\mathrm{log}\left(u+\frac{X_0}{v}\right)$$
(20)
with $`X_0`$ the position of the mirror at $`t=0`$. The corresponding normalised mode basis is given by

$`v_n={\displaystyle \frac{1}{2\sqrt{n\pi }}}\left[e^{i\mathrm{\Omega }_n\mathrm{log}(X_0/v+tx)}e^{i\mathrm{\Omega }_n\mathrm{log}(X_0/v+t+x)}\right],`$ (21)
$`\mathrm{\Omega }_n={\displaystyle \frac{n\pi }{\mathrm{arctanh}(v)}}n=1,\mathrm{},\mathrm{}`$ (22)
This set of exact analytical solutions proved useful in testing the accuracy of the numerical algorithm for the moving boundary problem, as detailed in Appendix A. An interesting feature of the solution is that for $`t_r=2X_0/(1v)`$ the evolved modes coincide with the original ones, re-scaled to the new box length. This feature is also present for the in-out motion studied above. It implies that at $`t_r`$ the in modes are identical to the out modes and the Bogolugov coefficients should vanish. As a consequence we should observe no particle production for this particular stopping time. We checked that this ‘resonance’ does indeed take place by performing the simulation described above with final time set to $`t_f=t_r`$. For several choices of expansion velocity (since $`t_r`$ grows with the inverse of $`1v`$ leading to very long simulation times, we restricted ourselves to the low-$`v`$ regime) the total number of particles produced was indeed zero.

### B Oscillating Mirror

Moving a step up in complexity from the uniformly moving mirror problem, we will consider in this section a cavity whose boundary executes small, periodic oscillations. This type of dynamical Casimir effect setting has been widely studied, and there is hope that an equivalent of a vibrating mirror could be experimentally realized. In the particular situation where the mirror frequency resonates with one of the cavity modes, particle production effects may be magnified, leading eventually to observable results. A thorough review of the topic with extensive references can be found in .

We will consider a sinusoidal mirror trajectory given by

$$X(t)=X_0+\frac{A}{2}\left[1\mathrm{cos}(\omega t)\right]$$
(23)
In the simulations described below, the amplitude of the motion $`A`$ will be taken to be considerably smaller than $`X_0`$, the initial size of the box. Typically we will set $`X_0=50`$ and $`A=2`$. The frequency of the motion $`\omega `$, will vary for different runs, care being taken that the velocity of the mirror never exceeds the speed of light. In general $`\omega `$ will be set to be one of the natural frequencies of the cavity, $`\omega _n=n\pi /X_0`$. This choice of frequencies is usually favoured in the literature, since it is expected that resonance between the cavity modes and the motion of the mirror will lead to higher particle production for long times. The resonant case is also considerably easier to tackle analytically (see  where the problem is solved exactly for $`n=2`$ and  for a calculation of $`R(u)`$ for a family of periodic-like trajectories). In fact, as far as we are aware, off-resonant solutions are only discussed in , where the energy production is evaluated for the general case. Our method of course, does not distinguish between resonant and off-resonant frequencies - in general we will present results for resonant frequencies and discuss how they change for the general case.

We start by looking at results for very short times, when the mirror executes a small number of oscillations. In Fig. 3 we show the total number of particles produced as a function of time, for several choices of the oscillation frequency. At every time-step, the numerical field profiles were contracted with the static basis for a stationary box of length $`X(t)`$. As discussed in the previous section this will lead to a convergent particle spectrum, but only for times when the mirror velocity is zero. For the sake of clarity, and with this caveat in mind, we still plot the total particle number for all time steps. For the number of modes simulated, the results for non-stationary final states do not deviate much from the finite case. Since the divergence is logarithmical, we presume it would only manifest itself if much higher frequencies were included.

It is clear from Fig. 3 that the evolution is effectively linear for the period shown. The production rate changes only for times of the order of twice the box length - up to then each oscillation produces roughly the same increase in particle number. This suggests that for small times, the field disturbances radiate away from the moving boundary without significant interference. For later times, after reflecting back on the stationary wall, they will interact with the vibrating mirror, constructively or otherwise depending on their frequency, leading to a change in the production rate.

For each value of the oscillating frequency $`\omega `$, we fitted the results to a linear law of the form $`N_{\mathrm{tot}}(t)=\gamma t`$. We checked that allowing for a constant term in the linear expression, $`N_{\mathrm{tot}}(t)=b+\gamma t`$ did not change the estimate for the production rate $`\gamma `$ ($`b`$ was always found to be very small, of order less than $`10^3`$). Clearly, higher oscillating frequencies lead to higher particle production rates. In Fig. 4 the dependence of $`\gamma `$ on the mirror’s frequency is made more explicit. The results are very well fitted by a power law of the form $`N_{\mathrm{tot}}(t)=D\omega ^kt`$ with $`k=3.3`$, a result in agreement with the analytical estimates in the literature. It is generally predicted that for a single oscillation the number of particles produced varies as a power law of the oscillation frequency. The particular value of the power depends on the specific type of oscillating trajectory. In  the production power for half an oscillation was shown to vary between $`2.`$ and $`4.`$, and in  a $`\omega ^4`$ dependence was found in the limit of small velocities. For a single oscillation, our result corresponds to a power of $`2.3`$, which is in line with these values.

As for the amplitude coefficient $`D`$, it was suggested in  that it should vary quadratically with the amplitude of the vibration of the mirror. We performed a quick test of this estimate by repeating the same set of runs with $`A`$, the mirror oscillating amplitude set to $`1.`$, half its previous value. Fitting both sets of data simultaneously to the same power $`k`$ and allowing $`D`$ to vary, we obtained $`k=3.2`$, $`D=0.0342`$ and $`D=0.00780`$ for $`A=2.`$ and $`A=1.`$ respectively. The rate of the two amplitudes is $`0.23`$, consistent with a reduction by a factor of $`4`$ as expected.

We now look at the longer time behaviour of the system. As mentioned above, the evolution stops being described by a linear law at about $`t2X_0`$, when a new regime with a different production rate sets in. In Fig. 5 we can see the total particle number for longer times for the resonant frequency $`\omega =2\pi /X_0`$.

The evolution for the period of time shown is well fitted by a power law $`N_{\mathrm{tot}}=8.5\times 10^7t^\alpha `$, with $`\alpha =1.85`$. For the dynamical range in question a power law growth with $`\alpha =2.`$ was predicted in ref. . This is in reasonable agreement with our numerical estimate.

For higher frequencies, as far as we are aware, there are no explicit predictions for the particle number time dependency (see  for a calculations of the energy produced). In Fig. 6 we show $`N_{\mathrm{tot}}`$ as a function of time for two higher resonant frequencies $`\omega =6\pi /X_0`$ and $`\omega =7\pi /X_0`$. We also include $`N_{\mathrm{tot}}`$ for the off-resonant case $`\omega =6.5\pi /X_0`$.

The overall evolution does not fit any simple expression in any of the three cases. We note that the initial period of linear evolution is always present and we have checked that the production rates follow the power law shown in Fig. 4, even in the off-resonant case. For later times though, the contrast between off-resonant and resonant trajectories is very clear, with limited quasi-periodic production in the former case and continuing increase in the later. It is interesting not to find any qualitative difference between the long time behaviour of the even and odd resonant frequencies systems, n=6 and n=7. In fact it has been often suggested that the main phenomena governing the evolution in the long time limit is akin to parametric resonance . If that were strictly the case though, resonant amplification would only take place for oscillating frequencies equal to twice the cavity’s natural frequencies . That is, using our notation, we should see long time suppression of the production rate for odd values of $`n`$, which clearly is not the case. This seems to confirm the results in  where it was found that the energy produced grows exponentially for long times, for any motion with $`\omega `$ equal to a multiple of the fundamental frequency $`\pi /X_0`$.

We now look at the distribution of the particles produced in terms of their frequency. In Fig. 7 we show the late time particle spectrum for the resonant $`n=6`$ motion. We observe an alternating succession of excited and suppressed frequencies. In particular we see that there is no production at all for frequencies multiple of the oscillating frequency. This cancellation is exact and it coincides strictly with points of the form $`N\omega `$, with $`N`$ integer. The peaks of the spectrum follow a similar but less precise pattern. The most excited frequency is $`\omega /2`$ and the following peaks take place at about points $`M\omega +\omega /2`$, though there seems to be a tendency for the peak to move up for higher frequencies. The fact that the maximum of the production rate corresponds to a frequency of about half the excitation, with higher even multiples suppressed and odd ones enhanced, is once again reminiscent of parametric resonance. Nevertheless, since the same pattern is observed for oscillating frequencies odd multiples of $`\pi /X_0`$, we must conclude that a more complex process is operating in the system.

Finally we contrast the resonant spectrum with the $`n=6.5`$ off-resonant case as shown in Fig. 8. Clearly, the overall amplitude is greatly reduced (by about one order of magnitude) and the pattern of alternating peaks is no longer present. The single peak is once again centred around $`\omega /2`$ but its height oscillates widely in time, leading to the periodic fluctuations in total particle number shown in Fig. 6.

## IV Non-uniform time-varying Dielectric

In the examples discussed above, the varying geometry of the systems considered was the primary cause behind particle production from the initial vacuum state. In this Section we will look at a system with fixed dimensions, but whose bulk physical properties change both in space and in time. In particular we will consider wave propagation in a background with a space and time dependent dielectric “constant”. The system will be modelled by a straightforward generalisation of the wave equation

$$\frac{}{t}\left(ϵ(x,t)\frac{\varphi }{t}\right)=\frac{^2\varphi }{x^2},$$
(24)
where $`ϵ(x,t)`$ is the varying dielectric term. Other generalisations are possible, their form depending ultimately on the properties of the particular system one chooses to model. In , the change in $`ϵ`$ is caused by a moving half-infinite dielectric material, and compliance with Lorentz invariance leads to a wave equation with a cross derivative term. In so-called slow light experiments , based on electromagnetically-induced transparency, the optical properties of the medium are controlled by an external source. These systems can be described by an effective Lagrangian that includes a mass-like term for the field, leading to yet another choice of field theory. Here, for the sake of simplicity we follow  and use Eq. (24), though our approach can be easily extended to the cases mentioned above.

For two solutions of Eq. (24) the time invariant scalar product is given by

$$(\varphi _1,\varphi _2)=i_0^Lϵ(x,t)\left[\varphi _1_t\varphi _2^{}_t\varphi _1\varphi _2^{}\right]𝑑x$$
(25)
the fields being defined in a finite box of size $`L`$. We will look at cases where $`ϵ(x,t)`$ is spatially uniform both for the initial and final times, $`t_1`$ and $`t_2`$. The corresponding normalised mode basis are given by:

$`u_n(x,t)=\sqrt{{\displaystyle \frac{c_i}{n\pi }}}e^{i\omega _nt}\mathrm{sin}(k_nx),`$ (26)
$`\omega _n=c_ik_n,k_n={\displaystyle \frac{n\pi }{L}}n=1,\mathrm{},\mathrm{}`$ (27)
where $`c_i=1/\sqrt{ϵ_i}`$ with $`i=1,2`$ denotes the initial and final speed of light respectively.

In the specific scenario we will be simulating, the original medium with dielectric constant $`ϵ_1`$ is replaced throughout the box by a medium characterised by $`ϵ_2`$. Medium $`1`$ expands into medium $`2`$ from left to right with a given velocity $`v`$, as a one-dimensional analogue of an expanding bubble. At intermediate time steps, the two regions with different dielectric constants are separated by a “wall” that interpolates between $`ϵ_1`$ and $`ϵ_2`$. We model the wall profile as a sine

$$\chi (x)=\{\begin{array}{cc}ϵ_1,& \text{if }x>d/2\hfill \\ ϵ_+ϵ_{}\mathrm{sin}(\pi x/d),& \text{if }d/2<x<d/2\hfill \\ ϵ_2,& \text{if }x<d/2\hfill \end{array}$$
(28)
with $`ϵ_+=(ϵ_1+ϵ_2)/2`$, $`ϵ_{}=(ϵ_1ϵ_2)/2`$, $`d`$ being the wall thickness. Using the wall profile $`\chi (x)`$ we can easily define a time-dependent dielectric function modelling the setting described above, as $`ϵ(x,t)=\chi (xvt)`$ with $`t`$ varying between $`t_1=d/(2v)`$ and $`t_2=d/(2v)+L/v`$. We will be looking at cases where the final speed of light is less than the original one, corresponding to a vacuum box being filled with a denser medium. The initial dielectric constant will be set to $`ϵ_1=1`$ with $`ϵ_2>ϵ_1`$. In general we will constrain the wall expansion velocity to be less than the speed of light in both the initial and final medium, that is $`v<1/\sqrt{ϵ_2}<1/\sqrt{ϵ_1}`$

In Fig. 9 we illustrate the dependence of the number of particles produced on the expansion velocity. The wall thickness was set to $`d=5.`$, considerably smaller than the physical box size, $`L=50.`$ The initial dielectric constant was fixed to $`ϵ_1=1.`$ whereas its final value $`ϵ_2`$ was varied between $`1.2`$ and $`4.0`$. For each particular choice of $`ϵ_2`$ we simulated a series of expansion velocities ranging from $`v0.3`$ up to the limit $`v1/\sqrt{ϵ_2}`$.

We found that for low velocities the production rate was considerably suppressed. As the expansion velocity increases, the total number of particles produced goes up steeply - on average by more than two orders of magnitude - as shown in Fig. 9. For mid and high values of the velocity, the dependence of the final particle number on $`v`$ is well described by an exponential law, $`N_{\mathrm{tot}}=Ae^{\gamma v}`$. This expression is valid up to values of $`v`$ considerably near the speed of light in the second medium.

When the final dielectric constant is allowed to vary we observe as expected, that stronger changes in the optical properties of the system lead to higher production rates. As a consequence we should observe an increase in the coefficient $`\gamma `$ for higher values of $`ϵ_2`$. This is fact the case, and it turns out that the dependence of $`\gamma `$ on the final dielectric constant is very simple. In the regime studied $`\gamma `$ varies linearly with $`ϵ_2`$ as can be seen in the fit shown in the right plot of Fig. 9. Using dimensionality arguments we can then write

$$N_{\mathrm{tot}}\mathrm{exp}\left[\left(K_0+K_1\frac{ϵ_2}{ϵ_1}\right)\frac{v}{c_1}\right]$$
(29)
where the dimensionless constants in the exponential are given by $`K_0=12.0`$ and $`K_1=5.1`$ for the system in question.

We also tested how the properties of the system depend on the thickness of the wall separating the two media. For a fixed value of the final dielectric constant $`ϵ_2=2.`$ we allowed $`d`$ to vary. A sample of the results is shown in Fig. 10. As expected, for larger values of $`d`$ the transition between the initial and final media is smoother, and the production rate is lower. As before, for $`d8`$ the final particle number depends exponentially on the expansion velocity, though for higher values of the wall thickness we observe deviations from this simple law. In the “thin wall” regime, the exponent $`\gamma `$ varies very little with $`d`$, as can be seen in the right plot in Fig. 10. For $`1d7`$ the change in $`\gamma `$ is less than 10%.

Finally we checked whether these results depend on the specific shape of the wall profile. We replaced the sin wall by one where a linear function interpolates between the two dielectric constants. We observed that the results were virtually unchanged.

## V Conclusions

We have developed and tested the numerical techniques necessary to tackle the problem of particle production from a vacuum state in general scenarios. In particular, we have shown that the approach performs quite well when applied to cosmology, to the dynamical Casimir effect and to systems with complex optical properties. In some of these cases we were able to use the numerical data to obtain simple phenomenological laws. This illustrates the power of the method and shows how it can be used as a tool to approach problems for which there is little or no analytical information. The opportunities are many, as was mentioned in the bulk of the paper. In the context of the dynamical Casimir effect, the first obvious step will be to generalise this approach to higher dimensions, where little is known analytically. This poses no obvious numerical difficulty, in contrast to the analytical problems raised by the loss of conformal invariance. Extensions to general types of geometry and the inclusion of a mass term in the theory should also be within reach of the method. Optical systems in $`2D`$ and $`3D`$ with alternative Lagrangians, and realistic phase transitions can also be subject to the same techniques. Finally, we should stress that since in general we are able to compute the full set of Bogolugov coefficients, we do not have to restrict the initial state to be in the vacuum. Time evolution of multi-particle states such as thermal configurations, can be easily simulated within this framework.

###### Acknowledgements.

N. D. A. was supported by PPARC. N. D. A. would like to thank Claudia Eberlein and Fernando Lombardo for useful suggestions. Part of the research was done within the framework of the E.S.F. COSLAB network.

## A Numerical Algorithm for Moving Boundary Conditions

In all simulations discussed above the scalar field wave equation was solved using a straight-forward Leap Frog algorithm (see for example ). Only the case of the moving mirror (Section III) requires more attention. For a box with time-dependent length, if the lattice spacing is kept constant, the number of points in the simulation domain must change and this may lead to large numerical errors.

The Leap Frog (LF) discretization of the wave equation Eq. (15) is

$`\mathrm{\Pi }_{i,n+1/2}=\mathrm{\Pi }_{i,n1/2}+dt^2\varphi _{i,n}`$ (A1)
$`\varphi _{i,n+1}=\varphi _{i,n}+dt\mathrm{\Pi }_{i,n+1/2},`$ (A2)
where $`i`$ is the spatial lattice index and $`n`$ the time step. As usual with LF-type algorithms, the field and its momentum are defined respectively on integer and half-integer time-steps. As before the coordinate of the moving boundary will be denoted by $`X(t)`$, and at each time-step its value will not necessarily lie exactly on one of the lattice spatial points. This discrepancy will have to be taken into account. In particular, as $`X(t)`$ increases new points will be introduced in the lattice and one must find a rule for assigning a field or momentum value to these.

Let us start with the case where the new lattice point $`i+1`$ lies in a half-integer time slice $`n+1/2`$, as depicted in Fig. 11.

We start by differentiating in respect to time the boundary condition $`\varphi (X(t),t)=0`$,

$$\frac{\varphi }{x}_{(X(t),t)}\dot{X}(t)+\mathrm{\Pi }(X(t),t)=0$$
(A3)
This allows us to calculate the momentum at the boundary in terms of the spatial derivative of the field. Since the field vanishes at the boundary, we can approximate $`_x\varphi `$ at the point $`(i+1,n+1/2)`$ by $`(\varphi _{i,n}+\varphi _{i,n+1})/(2dx)`$, using the fact that $`\varphi _{i,n+1}`$ can be evaluated in advance. This leads to:

$$\mathrm{\Pi }_{i+1,n+1/2}=\frac{\dot{X}_{n+1/2}}{2dx}(\varphi _{i,n+1}+\varphi _{i,n})$$
(A4)
Note however that here we implicitly assume that the boundary crosses the lattice at precisely the point $`(i+1,n+1/2)`$, which is not necessarily true. The above estimate can be improved by introducing an auxiliary fictitious lattice point $`(i+1,n+1/2\beta )`$, the exact location where the mirror’s trajectory intersects the lattice. The time difference between the crossing and time-step $`n+1/2`$ is given by $`\beta dt`$ and can be evaluated by interpolating the boundary trajectory between steps $`n`$ and $`n+1/2`$. At the crossing point we have

$$\mathrm{\Pi }_{i+1,n+1/2\beta }=\frac{\dot{X}_{n+1/2}}{2dx}\left[(12\beta )\varphi _{i,n+1}+(1+2\beta )\varphi _{i,n}\right]$$
(A5)
Note that by setting $`\beta =0`$ we recover the previous estimate Eq. (A4). Using Eq. (A5) we can then update the field at time-step $`n+1`$ as

$$\varphi _{i+1,n+1}=dt(1/2+\beta )\mathrm{\Pi }_{i+1,n+1/2\beta }$$
(A6)
The next momentum up-date, at $`n+3/2`$, will be calculated taking into account that boundary was placed at $`(i+1,n+1/2\beta )`$.

This method is easily generalised to the situation where the newly introduced lattice point lies on a integer time slice $`n+1`$. We proceed in a similar way, defining an auxiliary momentum point at $`(i+1,n+1/2\beta )`$, though in this case $`\beta `$ will be negative. The following updates proceed exactly in the same way as described in Eqs. (A5) and (A6).

When the trajectory of the mirror is such that no new points need to be introduced, we must still take some care when updating the fields near the boundary. In particular, the calculation of the second derivative term $`^2\varphi _{i,n}`$ in Eq. (A2) should reflect the fact that the boundary may lie between lattice sites. Referring again to Fig. 11 we define $`\alpha `$ as the distance in lattice units from point $`(i,n)`$ to the mirror. Since the field vanishes at the boundary we will have

$$^2\varphi _{i,n}=\frac{\alpha \varphi _{i1,n}(1+\alpha )\varphi _{i,n}}{\alpha dx^2}$$
(A7)
By setting $`\alpha =1`$, this expression reduces to the less accurate estimate for the second derivative, where the boundary’s position is identified with lattice point $`i+1`$.

In order to check the accuracy of the algorithm, we have used the exact analytical solution Eq. (22) for the uniformly expanding box and compared it to the corresponding numerical solution. In Fig.12 we show for a particular solution, the numerical errors for three different methods based on the approximations discussed above. The error is defined as the (normalised) integral of the square of the difference between the exact and numerical solutions. Algorithm $`1`$ is obtained by setting $`\beta =0`$ and $`\alpha =1`$ in Eqs. (A5) and (A7). This corresponds to a straight forward implementation of the Leap Frog method with no smoothing at the boundary in either time or space directions. Method $`2`$ includes the corrections for the momentum at the boundary but non-improved spatial derivatives ($`\alpha =1`$). Finally, in method $`3`$ both $`\alpha `$ and $`\beta `$ are evaluated at each time-step according to the mirrors’s position. The improvement in accuracy is clear: algorithm $`3`$ decreases the errors by a factor of $`10`$ over method $`2`$ and by a factor of $`25`$ over method $`1`$. For the purposes of the simulations in this work this is enough. Direct observation of the numerical solutions shows that the moving boundary introduces high-frequency discontinuities in the fields (in particular in the momentum). These are considerably reduced by method 3, and in the end they tend to average out when the Bogolugov coefficients are evaluated, leading to reasonably small errors in the final physical quantities.