format: text

text:
# Goldstone excitations from spinodal instability

## 1 Introduction

This investigation is the first stage of a complex study with the aim to explore in detail the mechanism of the energy transfer to the Higgs field at the end of the inflationary period of the Universe as it is anticipated in scenarios of hybrid inflation .

The $`O(2)`$-symmetric model used in studies of hybrid inflation has the following Lagrangian density for the real scalar (“inflaton”) field $`\psi `$ and the complex (“Higgs”) field $`\mathrm{\Phi }=\mathrm{\Phi }_1+i\mathrm{\Phi }_2`$:

$$L=\frac{1}{2}(_\nu \psi (x))^2\frac{1}{2}m_\psi ^2\psi ^2+\frac{1}{2}|_\nu \mathrm{\Phi }(x)|^2\frac{1}{2}m_\mathrm{\Phi }^2|\mathrm{\Phi }|^2\frac{\lambda }{24}|\mathrm{\Phi }|^4\frac{1}{2}g^2\psi ^2|\mathrm{\Phi }|^2,$$
(1)
where $`m_\mathrm{\Phi }^2`$ has the wrong sign.

At the beginning the energy density is carried partly by the inflaton field due to its $`𝒪(M_P)`$ initial amplitude and partly by the potential energy of the Higgs-field which starts at the symmetric maximum. It is transformed very efficiently into kinetic and gradient energy densities when the squared effective mass of the Higgs field becomes negative and the modes of the complex field with low spatial frequency $`|𝐤|<|M_{eff}|`$ start to increase exponentially due to the spinodal instability . The Higgs field eventually arrives to a symmetry breaking ground state, on its top with massive and massless thermal excitations.

The present study was carried out in Minkowski metrics and with classical fields starting from an initial state corresponding to the above situation. Therefore at this stage the results have more relevance to non-equilibrium phase transitions in relativistic field theoretical models than to cosmology. The results will serve as reference for future simulations to be performed in an expanding FRW-geometry. Our study was inspired by recent papers of Felder et al. describing also the evolution of a complex Higgs field after a sudden change of sign of its squared mass .

Very recently a detailed study concentrating on the statistics of defect formation was published , where the change of sign in the mass of the Higgs field is induced by the free motion of a homogenous real scalar field. The present work represents a complementary approach since we study rather the elementary excitations of the system than extended objects.

Our aim is to investigate systematically the way the different components of the complex Higgs field are excited. Special attention is paid to the real time appearance of the Goldstone mode after the period of spinodal instabilities. Analogous question was analyzed recently by Boyanovsky et al.  and by Baacke and Heitmann  in the infinite component (large $`N`$) limit of the quantum dynamics. They solved the coupled set of the renormalised equations of motion of the order parameter and of the to-be-Goldstone modes following a squared mass quench. The effective mass square governing the mode equations of the Goldstone fluctuations was analyzed which depends non-linearly on the actual value of the order parameter and reflects also the global back-reaction of the Goldstone-modes. It was found that the squared mass of the ”pion” modes starts to oscillate around the asymptotic zero value rather early, but the oscillations are damped only for asymptotically large times as $`t^1`$. In a recent paper  renormalised non-equilibrium gap equations were solved for the effective mass-squares of the longitudinal and transversal modes propagating on a time dependent background. These masses were in earlier papers  with a self-consistent parametrisation of the corresponding propagators. It turned out that the time-dependent variational mass-squares do not obey the Goldstone theorem away from equilibrium. It will be interesting to compare these semi-analytical approximate investigations with the direct mass measurements performed in numerical simulations of $`O(N)`$ symmetric systems.

The classical (cut-off) field theory provides a useful point of reference, since gapless excitations are present in its broken symmetry phase near equilibrium. Our numerical results for $`N=2`$ hint at an essentially different picture on the real time genesis of the Goldstone modes.

The presentation of our results is organised in the following way. Section 2 summarizes the set of parameters used in and the algorithmical details applied to the numerical solution. Our detailed discussion is divided into two parts. In section 3 the methods for finding the independently moving degrees of freedom are presented. There we give direct evidence for the early presence of massless excitations. In section 4 the process of thermalisation is described, pointing out the slow damping of angular O(2) oscillations, supporting the presence of Goldstone modes. Section 5 contains our conclusions.

## 2 Lagrangian parameters, initial conditions and discretisation

The parameters we choose for the present investigation imitate a situation which would be characteristic for a GUT-like transition. In lattice units the Higgs particle has unit squared mass parameter (with the wrong sign), the inflaton mass square is $`10^3`$ times smaller. The initial amplitude of the inflaton is large: $`\psi _0=15`$. Variation of its value does not influence the qualitative features of the time evolution of the system until its contribution to the potential energy is negligible relative to the initial energy content of the Higgs field.

A very weak nonlinearity $`\lambda =24\times 10^4`$ was fixed in order to facilitate the comparison of our results with . In this way the dominant part of the initial energy was stored in the Higgs potential. When we varied $`\psi _0`$ in a range, without changing the dominant component of the potential energy, no change was experienced in the thermalisation scenario after the spinodal instability. The use of a stronger Higgs self-coupling does not seem to have any qualitative consequence on the features of the appearance of Goldstone modes.

The initial space average of the Higgs field was fixed to zero. To the two real components of the complex field at every site a (white noise) amplitude evenly distributed on the $`[1,1]\times 10^4`$ interval was assigned. The amplitude of this random noise can be used to control which phase the system finally relaxes to. With the present choice the final value of $`|\mathrm{\Phi }|`$ is very close to the classical minimum of the Higgs potential, the system is very cold. One needs 5 order of magnitude stronger noise to drive the system into the symmetric phase.

Runs starting with different initial noise seeds were used to assess the statistical errors of the mass-values extracted. The restricted number of runs (6-14 depending on the size of the system), might cause some underestimation of the errors.

The value of the inflaton-Higgs coupling used in the present investigation is $`g^2=10^2`$. It turned out in our numerical studies that with stronger coupling (for $`g^28\times 10^2`$) the time necessary for the full development of the spinodal instability takes several oscillation periods of the inflaton, and the growth of unstable modes can happen only during the passage of the inflaton field through zero. This behaviour can be understood with help of the theory of broad parametric resonance developed in  when it is applied to a system where also spinodal instability can occur. The main reason for this phenomenon is the fact that with the increase of $`g^2`$ the time interval of the non-adiabatic inflaton evolution becomes shorter. Its discussion is beyond the scope of the present paper. It has been checked that the spectral characterisation of the different modes as they appear after the spinodal instability are independent of $`g^2`$.

We have followed the temporal evolution of the inflaton as well as of the $`O(2)`$-field on a spatial lattice with the equations written for dimensionless field quantities:

$`\ddot{\psi }(x)\psi (x)+{\displaystyle \frac{m_\psi ^2}{m_d^2}}\psi (x)+g^2\mathrm{\Phi }^2\psi =0`$
$`\ddot{\mathrm{\Phi }}_i(x)\mathrm{\Phi }_i(x)+{\displaystyle \frac{m_\mathrm{\Phi }^2}{m_d^2}}\mathrm{\Phi }_i(x)+{\displaystyle \frac{\lambda }{6}}|\mathrm{\Phi }(x)|^2\mathrm{\Phi }_i(x)+g^2\psi ^2\mathrm{\Phi }_i=0`$ (2)
with the scaled quantities

$`{\displaystyle \frac{m_\mathrm{\Phi }^2}{m_d^2}}=1,{\displaystyle \frac{m_\psi ^2}{m_d^2}}=0.001,`$
$`{\displaystyle \frac{\psi _d}{m_d}}=\psi ,{\displaystyle \frac{\mathrm{\Phi }_d}{m_d}}=\mathrm{\Phi },t_dm_d=t,a_dm_d=a1`$ (3)
(the symbols with index $`d`$ refer to dimensionfull quantities).

The characteristic time scales related to the elementary excitations turned out much shorter than the relaxation times, especially for the long wavelength Goldstone modes. This circumstance made the choice of the time-step in the standard leapfrog algorithm used by us a sensitive issue. The stability of the algorithm was ensured by choosing the time step in the interval $`0.125a_t/a_s0.0625`$. The long term stability of the time evolution was controlled by monitoring the conservation of the energy and the $`O(2)`$ charge of the system. Due to the initial conditions the latter was fixed at zero. The charge density for the above time step was fluctuating around zero with an amplitude $`0.04`$, but in selected runs using hundred times shorter time steps we could reach an amplitude $`𝒪(10^4)`$ with unchanged qualitative features.

## 3 Finding the independent excitations

The spinodal instability driving the Higgs field from the symmetric phase to the broken symmetry state was investigated on lattices with increasing $`L^3=32^3,64^3,128^3`$ volumes and with periodic spatial boundary conditions. Due to our initial “white noise” amplitude distribution the system chooses a well-defined (though random) direction in the $`O(2)`$ internal space contrary to the case studied by  where the space average of the Higgs field stayed at the origin.

The coincidence of the change of sign of the instant effective squared mass $`M_{eff}^2`$ with the start of the spinodal instability is investigated in Fig.1. Various alternative definitions of the effective Higgs mass are used:

$`M_{eff,1}^2`$ $`=`$ $`m_\mathrm{\Phi }^2+g^2\psi ^2(t),`$
$`M_{eff,2}^2`$ $`=`$ $`M_{eff,1}^2+g^2(\psi \psi )^2,`$
$`M_{eff,3}^2`$ $`=`$ $`M_{eff,2}^2+{\displaystyle \frac{\lambda }{2}}|\mathrm{\Phi }|^2.`$ (4)
The brackets mean spatial averaging.

The first definition, where one assumes that the fluctuations of the $`\psi `$-field are negligible is used in most discussions of hybrid inflation. It predicts a change of sign when the inflaton reaches the critical value $`\psi _c=|m_\mathrm{\Phi }|/g`$, and was ensured to be valid in  by setting $`m_\psi ^2=0`$. In the present case the inflaton rolls in a harmonic potential and one can test if the fluctuations of the $`\psi `$ and $`\mathrm{\Phi }`$ fields have any effect at the start of the spinodal instability. In the Figure only $`M_{eff,1}^2`$ and $`M_{eff,3}^2`$ are displayed, since $`M_{eff,2}^2`$ is indistinguishable from the first one. The deviation of the two versions of the mass appears near the end of the instability, when $`M_{eff,1}^2`$ stabilizes at the value $`1`$, while $`M_{eff,3}^2`$ approaches a positive stationary value after wild oscillations. The increase of the order parameter $`\rho d^3x|\mathrm{\Phi }|/V`$ becomes noticeable with a slight time delay after the sign-change of the effective Higgs mass. This is due to the gradual increase of $`M_{eff,1}^2`$. The spinodal increase of the order parameter is exponential, with a slowly varying slope. Near the saturation of the order parameter the slope is very close to the unit value expected with our choice of $`m_\mathrm{\Phi }^2`$.

Next, we turn to the separation of the independently oscillating degrees of freedom. At every lattice point the complex field can be represented as $`\mathrm{\Phi }(𝐱,t)=\rho (𝐱,t)\mathrm{exp}(i\phi (𝐱,t))`$. The main question of this investigation was, at what time the radial and angular Higgs variables become dinamically independent. This can be monitored by calculating the temporal variation of the elements of the “velocity”-correlation matrix $`W_{ij}(t)`$, which was constructed the following way.

Let us define the 3-component quantity: $`v_i(𝐱,t)=\{\dot{\rho }(𝐱,t),\rho (𝐱,t)\dot{\phi }(𝐱,t),\dot{\psi }(𝐱,t)\}`$ and form the spatial fluctuation matrix

$$W_{ij}=\frac{\stackrel{~}{W}_{ij}}{\stackrel{~}{W}_{ii}^{1/2}\stackrel{~}{W}_{jj}^{1/2}},\stackrel{~}{W}_{ij}(t)=\frac{1}{V}d^3xv_i(𝐱,t)v_j(𝐱,t)\frac{1}{V}d^3xv_i(𝐱,t)\frac{1}{V}d^3xv_j(𝐱,t).$$
(5)
In Fig.2 we show the early time variation of the $`W_{12}(t)`$ element on three different size lattices. One observes that the level of Higgs-Goldstone cross-correlation decreases by about two orders of magnitude from lattice size $`L=32`$ to $`L=128`$. All non-diagonal correlation coefficients show maxima during the spinodal instability. In case of $`W_{13}(t)`$ it was found that the maximum of the inflaton-Higgs cross-correlation coefficient is independent of the size and is $`𝒪(1)`$. In case of $`W_{12}`$ and $`W_{13}`$ also the maximum decreases with the size of the system. It is unclear what is the thermodynamical limit of this maximum reached during the spinodal instability.

After the saturation of the unstable modes the average level of all off-diagonal elements of the normalized correlation matrix stabilizes for $`L=128`$ at $`𝒪(10^7)`$. This level of decoupling is reached the slowest for $`W_{12}`$. Therefore the radial and angular motion of the $`O(2)`$-field in the internal space can be considered with very good accuracy as independent starting from about $`t200`$. The analysis gives very similar results when applied directly to the polar coordinates: $`x_i=\{\rho ,\phi ,\psi \}`$.

This analysis provides the basis for a separate application of various methods of mass determination to each of the three independent coordinates. Would the cross-correlations be $`𝒪(1)`$ beyond the instability, one could not avoid the application of a complicated diagonalisation process before applying the methods of dispersive mass determination to be described next.

The method developed and tested by us on the example of the one-component scalar field  has been further improved for the non-trivial task of the investigation of Goldstone modes. Its present algorithm is summarised in the following.

One keeps track of the time evolution of the lower-$`k`$ spatial Fourier transform of all three $`v`$-type lattice fields. No binning is applied. Next, one computes the frequency-spectra of the time evolution of the power of the single $`k`$-modes. They show several peaks and the strongest one was identified with the eigenfrequency of the investigated $`𝐤`$-mode (the others correspond to various forced oscillation modes). This correspondence yields the dispersion relation $`\omega (𝐤)`$.

In equilibrium systems we verified the correctness of this method by taking long time intervals. In Fig.3 such spectra are displayed for the lowest two $`𝐤`$ modes for both the Higgs and the Goldstone variable on the largest $`L=128`$ lattice. We observed essentially no volume dependence. The straight line extrapolation of the dispersion relation from finite $`𝐤`$ modes towards zero wave number is of very good quality and shows the absence of any gap in case of the angular mode. The Higgs-mode follows a massive relativistic dispersion relation.

It was checked that the analysis of shorter time intervals down to $`T=32|m_\mathrm{\Phi }|^1`$ results in masses which are in good agreement with spectra obtained with much finer $`\omega `$ resolution.

Far from equilibrium fields vary fairly fast, any analysis averaged over a long time interval misses the changes which occur on short time scales. With the present choice of the parameters a fast Fourier analysis (FFT) using $`T/|m_\mathrm{\Phi }|=32`$ allows an accurate enough determination of the Fourier coefficients in the interval $`\omega (2\pi /T,12\pi /T)`$. 14 independent runs were analysed on an $`L=128`$ lattice. The masses were obtained by subtracting the lattice $`|𝐤|^2`$ from the corresponding $`\omega ^2`$ for 9 modes, and fitting the remaining numbers to a constant. The average standard deviation of the fits performed in different moments is shown as an error estimate of the mass determination. Fig.4 displays the time evolution of the independent inflaton, Higgs and Goldstone masses.

In an earlier version of this method  we have assumed that all $`𝐤`$-modes oscillate with a single eigenfrequency. Then from the ratio $`|\dot{x}_i(𝐤)|^2/|x_i(𝐤)|^2`$ formed with help of $`x`$-type variables one can deduce the dominant frequency. In order to minimize fluctuations the numerator and the denominator were averaged separately over the modes lying in narrow $`|𝐤|`$ bins. The squared frequency was linearly extrapolated to $`|𝐤|=0`$ and the limiting value identified with the squared mass.

In Fig.5 one can see that the masses determined with this method provide for both components slightly higher values than the improved method. Their temporal fluctuations quickly diminish with increasing lattice size. The mass estimates shown in these figures represent averages of 40 independent runs for $`L=32`$ (9 runs were analysed for $`L=64`$, and 6 for $`L=128`$ with no sizable finite size dependence). The error bars displayed at the right hand side of both figures represents the standard deviation of the squared masses. Within these limitations on the accuracy the results of the two mass determinations are compatible.

For the Goldstone mode the best fit was reached by analyzing the phase factor $`|\dot{X}_k|^2/|X_k|^2,X=\mathrm{exp}(i\phi )`$. This choice was suggested by the idea of non-linear $`\sigma `$-models. We note that only modes with $`|𝐤|0.5`$ were included into the analysis. As we will discuss in the next section for lower wave numbers the Goldstone modes are highly excited during and after the instability. They relax very slowly, and their motion does not seem to follow a single eigenfrequency. From the right hand picture of Fig.5 one recognizes an apparent (rather small) negative squared mass fit for a short time interval directly after the spinodal instability. For later times the fit is compatible with vanishing mass.

We conclude this section with the statement that Goldstone modes are present in the excitation spectra of the model very early. The important question to be discussed in the next section is the degree of excitation of this mode relative to the others, that is the part of energy carried by the light modes.

## 4 Thermalisation and Goldstone-damping

The temporal evolution of the kinetic power spectra reflects very well the approach to classical thermal equilibrium, characterized by the equipartition of the energy. In Fig.6 the flattening of the kinetic Goldstone power spectra in the region $`k1`$ is clearly observed. The low $`𝐤`$ peak has an almost time independent height in the interval $`t=(100,1500)`$, its decrease is well recognizable only for $`t10^5`$.

This observation is probably the most appealing result of the present investigation. Below we attempt its semiquantitative interpretation with help of a combination of perturbative and non-perturbative arguments.

In a recent publication (, Appendix B and C) we have presented the perturbative evaluation of the Goldstone and Higgs damping rates in the broken phase of a classical(!) $`O(N)`$ symmetric field theory. The analysis was performed for constant value of the order parameter. The Goldstone and Higgs modes were characterized by constant mass values and by uncorrelated spectral functions

$$\mathrm{\Delta }_a(p)=2\pi ϵ(p_0)\delta (p^2M_a^2)f_a(p_0),a=H,G.$$
(6)
The induced source densities for these equations were found iteratively, assuming weak non-linearities in the system. The “self-energy” terms of the two types of excitations were established by taking into account the solution with just one iteration. The following explicit expressions were found:

$`\mathrm{\Pi }_H(p)`$ $`=`$ $`(\lambda \mathrm{\Phi })^2\left[R_{HH}(p)+{\displaystyle \frac{N1}{9}}R_{GG}(p)\right],\mathrm{\Pi }_G(p)=\left({\displaystyle \frac{\lambda \mathrm{\Phi }}{3}}\right)^22R_{GH}(p),`$
$`R_{ab}`$ $`=`$ $`{\displaystyle \frac{d^4q}{(2\pi )^4}\left(\mathrm{\Delta }_b(qp/2)\mathrm{\Delta }_a(q+p/2)\right)\frac{1}{2qpM_b^2+M_a^2}}.`$ (7)
In equilibrium systems the classical “number”-distribution is given by $`f_a(q_0)=T/q_0`$. The damping rates are determined by the imaginary parts of the self-energies.

The initial conditions generated by the spinodal instability for the later evolution is hard to specify explicitly. In principle one can attempt the numerical determination of the spectral function of each of the independent degrees of freedom . Here we followed a less systematic path based on the phenomenological description of the actual system presented in the preceding section, which confirms the presence of separate Higgs and Goldstone degrees of freedom with well-defined masses in our $`O(2)`$ invariant system shortly after the spinodal instability is over. Therefore we accept for the spectral densities the form of a free gas for each of them. Then the only remaining task is to find the nonequilibrium generalisation of the distribution functions $`f_a`$.

A classical analogue of the average occupation number was extensively used for the dynamical coordinate $`X_a`$ in :

$$_a(\omega _k)=\frac{1}{\omega _k}(|\dot{X}_a(𝐤)|^2+\omega _a(𝐤)^2|X_a(k)|^2),\omega _a(k)^2=𝐤^2+M_a^2.$$
(8)
It was computed here (see Fig.7) for the absolute value and the phase factor of the complex field. We see that the “number density” of the Higgs degree of freedom which initially decreases exponentially at later times approaches uniformly the $`1/\omega _k`$ regime. Also the virial equilibrium of the corresponding kinetic and potential contributions is fulfilled very early (for $`t500`$). In case of the Goldstone degree of freedom the convergence to the classical “number” density is not uniform, the excitation of the low-$`k`$ modes remains strong, and the approach to the classical regime in the high-$`k`$ region is slower.

The smoothness of $`_a`$ and the thermalisation tendency observed in their evolution leads us to the conjecture that the form (6) of $`\mathrm{\Delta }_a(p)`$ remains valid in the present case, only the functions $`f_a`$ are replaced by a far-from-equilibrium form. With this assumption we can use the method of evaluation described in Appendix C of . (One should be aware of the fact that beyond the present one-loop estimation the use of the proposed non-equilibrium form of $`\mathrm{\Delta }_a(p)`$ would lead to pinch singularities ). In case of the Goldstone modes we find the following expression for the on-shell damping rate:

$$\mathrm{\Gamma }_G(k_0=|𝐤|)=\left(\frac{\lambda \mathrm{\Phi }}{3}\right)^2\frac{1}{32\pi |𝐤|^2}_{M_H^2/4|𝐤|}^{\mathrm{}}𝑑p_0[_G(p_0)_H(p_0+k_0)].$$
(9)
The lower limit is a direct kinematical consequence of the zero mass nature of the Goldstone mode. It restricts the contribution to the large frequency tails of both distributions if the damping rate of the low-$`k`$ modes is to be estimated. Since these modes are filled very late after the spinodal instability is complete, the damping of these Goldstone modes is very inefficient.

It is amusing that the far-equilibrium classical evolution under specific initial conditions produces a phenomenon for the Goldstone excitation, which occurs at equilibrium in the quantum version of the model due to the exponentially small high frequency tail of the Bose-Einstein factor .

Due to the rather strong initial excitation of the low-$`k`$ Goldstone modes, the average kinetic energy stored in this degree of freedom is the biggest directly after spinodal instability. This averaging, however, is somewhat misleading. By subtracting the energy of the modes participating in the low-$`k`$ peak, the temperature of the ”equlibrated” modes appears somewhat lower than the corresponding Higgs-temperature. With time a metastable ”equilibration” between these modes becomes almost complete, as one can see in Fig.8.

This means that the equilibration in the $`O(2)`$ configuration space is much more efficient, than the energy transfer to the inflaton. Apparently, the inflaton field is frozen at a much lower temperature than the $`O(2)`$ fields.

The time spent by the inflaton in frozen state is of the same order of the magnitude as the relaxation time of the Goldstone low-$`k`$ peak. We have checked that with the increase of $`g^2`$ the inflaton temperature starts to increase earlier, but the transition becomes smoother and the equilibrium is reached later. On the other hand when the complex field is replaced by a real one ($`g^2`$ and $`\lambda `$ are kept at the same values) the equilibration time is an order of magnitude shorter. It is an interesting question if the final relaxation of the complex field leading to an equilibration with the inflaton was started by some sort of instability or it can be characterised by a smooth (exponential or power law) relaxation.

## 5 Conclusions

In this paper we have shown, that a well-defined Goldstone degree of freedom appears in an $`O(2)`$ invariant system, almost instantly after the spinodal instability, driven by the evolution of a real field coupled to it, is saturated. The thermalisation takes anomalously long due to the very efficient spinodal excitation of these modes which showed somewhat unexpectedly extremely slow relaxation in the low $`k`$ region.

This feature was understood with help of a leading order iterative estimation of the Goldstone damping rate relying on the existence of a far from equilibrium initial spectral function. The fact that a computation making explicit use of a heavy and of a massless field could explain correctly this phenomenon presents further evidence for the very early presence of the gapless modes which is the main result of this investigation.

The other surprising effect is the long (metastable) freezing of the real scalar field $`\psi `$ coupled to the Higgs mode relatively strongly. We provided numerical evidence that this phenomenon is also related to the presence of the Goldstone modes, but further study is needed for a complete understanding of the underlying mechanism.

It will be important to repeat the complete study in an expanding space-time geometry. Recent developments along this line include the computation of scalar particle production  and the dynamics of coupled real scalar fields  in an expanding Universe. The efficiency of the late equilibration of the inflaton is expected to depend critically on the Hubble parameter of the evolution.

During the spinodal instability the low $`k`$ modes become highly excited, the classical description is well justifiable for them. However, the transfer of the energy towards the stable high frequency modes might be more sensitive to quantum effects. Therefore in a second stage of the investigation above the upper cut-off $`K`$ applied to the modes treated classically $`(K>|m_\mathrm{\Phi }|)`$, quantum modes will be introduced and treated with help of renormalised mode equations .

## Acknowledgements

The authors gratefully acknowledge a valuable suggestion of D.J. Schwarz and careful reading of the manuscript by A. Jakovác and Zs. Szép. This investigation was supported by the Hungarian Research Fund (OTKA).


text:
# Experimental Implementation of Quantum Random Walk Algorithm

## I Introduction

Since the discovery of the first two quantum algorithms, Shor’s factoring algorithm 1  and Grover’s database search algorithm 2 , research in the new born field of quantum computation exploded 3 . Mainly motivated by the idea that a computational device based on quantum mechanics might be (possibly exponentially) more powerful than a classical one 4 , great effort has been done to investigate new quantum algorithms and, more importantly, to experimentally construct a universal quantum computer. However, finding quantum algorithms is a difficult task. It has also been proved extremely difficult to experimentally construct a quantum computer, as well as to carry out quantum computations.

Recently, several groups have investigated quantum analogues of random walk algorithms 5 ; 6 ; 7 ; 8 ; 9 ; 10 , as a possible direction of research to adapt known classical algorithms to the quantum case. Random walks on graphs play an essential role in various fields of natural science, ranging from astronomy, solid-state physics, polymer chemistry, and biology to mathematics and computer science 11 . Current investigations show that quantum random walks have remarkably different features to the classical counterparts 5 ; 6 ; 7 ; 8 ; 9 ; 10 . The hope is that a quantum version of the random walk might lead to applications unavailable classically, and to construct new quantum algorithms. Indeed, the first quantum algorithms based on quantum walks with remarkable speedup have been reported al . Further, experimental schemes have also been proposed to implement such quantum random walks by using an ion trap quantum computer 9  and by using neutral atoms trapped in optical lattices 10 . Up to now, only three methods have been used to demonstrate quantum logical gates: trapped ions 12 , cavity QED 13  and NMR 14 . Of these methods, NMR has been the most successful with realizations of quantum teleportation 15 , quantum error correction 16 , quantum simulation 17 , quantum algorithms 18 , quantum games 19  and others 20 . In this paper, we present the experimental implementation of the quantum random walk algorithm on a two-qubit NMR quantum computer, and we believe it is the first experimental implementation of such quantum algorithms.

We consider continuous-time random walks (CTRW) 6  rather than discrete-time random walks 7 , on a circle with four nodes. We show that the evolution of this quantum walk is periodic and reversible, and yields an exactly uniform probability distribution at certain time. While the classical CTRW is irreversible and only approximates the uniform distribution at infinite-time limit. Further, we experimentally implement this quantum walk on a two-qubit quantum computer, using a unitary operator which has the right “effective” Hamiltonian, with good agreements between theory and experiment. It is interestingly found that the property of the quantum walk strongly depends on the entanglement between the two qubits. The uniform distribution could be obtained only when the two qubits are maximally entangled.

## II Quantum CTRW on a Circle

The concept of quantum CTRW is proposed in Ref. 6 . On a circle with four nodes, we denote the set of nodes by $`\{0,1,2,3\}`$. Since the structure of the circle is periodic, only two nodes, $`\left(k+1\right)\mathrm{mod}4`$ and $`\left(k1\right)\mathrm{mod}4`$, are connected to node $`k`$ $`\left(k=0,1,2,3\right)`$. In the classical CTRW, let $`\gamma `$ denote the jumping rate, which is a fixed, time-independent constant that represents the probability moving from a given node to a connected one per unit time. The generator matrix 6  of this walk can therefore be written as

$$H=\left(\begin{array}{cccc}2\gamma & \gamma & 0& \gamma \\ \gamma & 2\gamma & \gamma & 0\\ 0& \gamma & 2\gamma & \gamma \\ \gamma & 0& \gamma & 2\gamma \end{array}\right).$$
(1)
Consider a particle walking classically on this circle, let $`P_k^C\left(t\right)`$ denote the probability of being at node $`k`$ at time $`t`$ (the superscript $`C`$ denotes “Classical”), then we have

$$\frac{d}{dt}P_k^C\left(t\right)=_{l=0}^3H_{kl}P_l^C\left(t\right).$$
(2)
This equation conserves the probability in the sense that $`_{k=0}^3P_k^C\left(t\right)1`$.

A natural way to the quantum version of this CTRW is to construct a Hilbert space spanned by the four basis $`\{|0,|1,|2,|3\}`$, respectively corresponding to the four nodes $`\{0,1,2,3\}`$. The state of a particle walking quantum mechanically on this circle is then denoted by $`|\psi \left(t\right)`$, which is generally a superposition of the four basis, rather than classical mixing of the probabilities of being at the four nodes. The generator matrix in classical walks is now treated as the Hamiltonian of the quantum evolution 6 . Therefore the Schrödinger equation of the state $`|\psi \left(t\right)`$ is

$$\frac{d}{dt}k|\psi \left(t\right)=i_{l=0}^3k\left|H\right|ll|\psi \left(t\right).$$
(3)
If measuring at time $`t`$, we can obtain a certain probability distribution on the circle. The probability of being at node $`k`$ is $`P_k^Q\left(t\right)=\left|k|\psi \left(t\right)\right|^2`$ (the superscript $`Q`$ denotes “Quantum”), and the conservation of probability is obviously guaranteed by the normalization $`\psi \left(t\right)|\psi \left(t\right)1`$.

Let the walking particle start from node $`0`$, it is then easy to find the probability of being at node $`k`$ at any time $`t`$ in both the classical walks and the quantum walks (by solving equations (2) and (3), respectively). The detailed calculation yields that for the classical walks the probabilities are

$$\{\begin{array}{c}P_0^C\left(t\right)=\frac{1}{4}+\frac{1}{2}e^{2\gamma t}+\frac{1}{4}e^{4\gamma t}\hfill \\ P_1^C\left(t\right)=P_3^C\left(t\right)=\frac{1}{4}\frac{1}{4}e^{4\gamma t}\hfill \\ P_2^C\left(t\right)=\frac{1}{4}\frac{1}{2}e^{2\gamma t}+\frac{1}{4}e^{4\gamma t}\hfill \end{array}.$$
(4)
In the quantum case, the initial state of the particle is $`|\psi \left(0\right)=|0`$, from equation (3), we have

$`|\psi \left(t\right)`$ $`=e^{iHt}|0`$
$`=e^{2i\gamma t}\mathrm{cos}^2\gamma t|0e^{2i\gamma t}\mathrm{sin}{}_{}{}^{2}\gamma t|2`$
$`+{\displaystyle \frac{i}{2}}e^{2i\gamma t}\mathrm{sin}2\gamma t\left(|1+|3\right).`$ (5)
Therefore the probabilities of the quantum walks are

$$\{\begin{array}{c}P_0^Q\left(t\right)=\mathrm{cos}^4\gamma t\hfill \\ P_1^Q\left(t\right)=P_3^Q\left(t\right)=\frac{1}{4}\mathrm{sin}^22\gamma t\hfill \\ P_2^Q\left(t\right)=\mathrm{sin}{}_{}{}^{4}\gamma t\hfill \end{array}.$$
(6)
The probabilities in equations (4, 6) are plotted in FIG. 1 (for the classical walks) and FIG. 2 (for the quantum walks) as functions of time $`t`$. From FIG. 1 and FIG. 2, we can see striking differences between quantum and classical random walks. FIG. 2 shows that the evolution of the quantum CTRW on a circle with four node are essentially periodic with a period $`T=\pi /\gamma `$. The particle walking quantum mechanically on this circle will definitely go back to its original position, and the evolution is reversible. It is also interesting to see that at time $`t=\pi /2\gamma `$ the probability distribution converges to node $`2`$. These phenomena are due to the quantum interference effects, which allows probability amplitudes from different paths to cancel each other.

To measure how uniform a distribution is, an immediate way is to use the total variation distance between the given distribution and the uniform distribution. In our case, the classical and quantum total variation distance as functions of time $`t`$ are

$`\mathrm{\Delta }^C\left(t\right)`$ $`={\displaystyle \frac{1}{2}}{\displaystyle _{k=0}^3}\left|P_k^C\left(t\right)\frac{1}{4}\right|,`$ (7)
$`\mathrm{\Delta }^Q\left(t\right)`$ $`={\displaystyle \frac{1}{2}}{\displaystyle _{k=0}^3}\left|P_k^Q\left(t\right)\frac{1}{4}\right|.`$ (8)
FIG. 3 depicts the dependence of $`\mathrm{\Delta }^C\left(t\right)`$ and $`\mathrm{\Delta }^Q\left(t\right)`$ on time $`t`$. From FIG. 3, we can see that the classical version of this walking process approaches the uniform distribution exponentially as time lapses. In contrast, the quantum process exhibits an oscillating behavior. An intriguing property of this quantum random walk is that $`\mathrm{\Delta }^Q\left(n\pi /4\gamma \right)=0`$ if $`n`$ is odd, which means that the probability distribution is exactly uniform at time $`t=\pi /4\gamma `$ and its odd multiples. While the classical walk can never reach the exactly uniform distribution, only approximates it at infinite-time limit.

## III Experimental Implementation

For the quantum CTRW on a circle with four node, the Hilbert space is 4-demensional. So it is natural to implement the quantum walks on a two-qubit quantum computer. The direct correspondence is to map the basis {$`|0`$, $`|1`$, $`|2`$, $`|3\}`$ of the quantum CTRW into the four computational basis {$`|0|0`$, $`|0|1`$, $`|1|0`$, $`|1|1`$}. This mapping is in fact to rephrase the number of nodes by the binary number system. Therefore the Hamiltonian in equation (1) can be written as

$$H=2\gamma II\gamma \left(I\sigma _x+\sigma _x\sigma _x\right),$$
(9)
where $`I`$ and $`\sigma _x`$ are the identity operator and the Pauli operator of a single qubit. The evolution operator of the two-qubit system is

$`U\left(t\right)`$ $`=e^{iHt}`$
$`=e^{2i\gamma t}\mathrm{exp}\left[i\gamma t\left(\sigma _x\sigma _x\right)\right]\mathrm{exp}\left[i\gamma t\left(I\sigma _x\right)\right].`$ (10)
And the state of a particle performing this quantum CTRW is

$`|\psi \left(t\right)`$ $`=e^{2i\gamma t}\mathrm{cos}^2\gamma t|00e^{2i\gamma t}\mathrm{sin}{}_{}{}^{2}\gamma t|10`$
$`+{\displaystyle \frac{i}{2}}e^{2i\gamma t}\mathrm{sin}2\gamma t\left(|01+|11\right).`$ (11)
It is interesting to investigate the relations between the distribution of the implemented quantum CTRW and the entanglement of the two-qubit state $`|\psi \left(t\right)`$. The entanglement of the two-qubit state $`|\psi \left(t\right)`$ in equation (11) can be directly calculated by the Von. Neumann entropy

$$S\left(t\right)=\mathrm{cos}^2\gamma t\mathrm{log}_2\left(\mathrm{cos}^2\gamma t\right)\mathrm{sin}{}_{}{}^{2}\gamma t\mathrm{log}_2\left(\mathrm{sin}^2\gamma t\right).$$
(12)
The correlation between the quantum total variation distance $`\mathrm{\Delta }^Q\left(t\right)`$ and the entanglement $`S\left(t\right)`$ is illustrated in FIG. 4. From FIG. 4, we can see that if there is no entanglement between the two qubits ($`S=0`$), $`\mathrm{\Delta }^Q`$ is at its maximal $`\mathrm{\Delta }^Q=\frac{3}{4}`$, which corresponds converging at node $`0`$ (or node $`2`$). While if the two qubits are maximally entangled ($`S=1`$), $`\mathrm{\Delta }^Q=0`$, which happens to be the situation that the walking particle is uniformly distributed on the four nodes. Therefore, we can say that the quantum random walk algorithm is enhanced by the quantum entanglement involved.

The quantum CTRW is implemented using our two-qubit NMR quantum computer. This computer uses a $`0.5`$ milliliter, $`200`$ millimolar sample of Carbon-13 labeled chloroform (Cambridge Isotopes) in d$`_\text{6}`$ acetone. In a magnetic field, the two spin states of $`^\text{1}`$H and $`^{\text{13}}`$C nuclei in the molecular can be described as four nodes of two qubits, while radio frequency (RF) fields and spin–spin couple constant $`J`$ are used to implement quantum network of CTRW. Experimentally, we performe twelve separate sets of experiments with various selection of time $`t`$ which is distinguished by $`\gamma t=\frac{n\pi }{12}\left(n=\{0,1,2,\mathrm{},12\}\right)`$. In the following, we replace the jumping rate $`\gamma `$ with $`\pi J`$ $`(J=215Hz)`$. In each set, the full process of the quantum CTRW is executed. We describe this experimental process as follows.

Firstly, prepare effective pure state $`|\psi \left(0\right)`$: The initial state in NMR is thermally equilibrium state $`\rho _{eq}4I_z^1+I_z^2`$ rather than a true pure state $`|\psi \left(0\right)`$. However, it is possible to create an effective pure state, which behaves in an equivalent manner. This is implemented as

$$R_x^1\left(\pi /3\right)G_zR_x^1\left(\pi /4\right)\tau R_y^1\left(\pi /4\right)G_z,$$

to be read from left to right, radio-frequency pulses are indicated by $`R_{axis}^{spins}\left(angle\right)`$, and are applied to the spins in the superscript, along the axis in the subscript, by the angle in the brackets. For example, $`R_x^1\left(\pi /3\right)`$ denotes $`\pi /3`$ selective pulse that acts on the first qubit about $`\widehat{x}`$, and so forth. $`G_z`$ is the pulsed field gradient along the $`\widehat{z}`$ axis to annihilate transverse magnetizations, dashes are for readability only, and $`\tau `$ represents a time interval of $`1/\left(2J\right)`$. Therefore, after the state preparation, we obatin effective pure state $`\rho \left(0\right)I_z^1+I_z^2+2I_z^1I_z^2`$ from equilibrium state $`\rho _{eq}4I_z^1+I_z^2`$.

Secondly, perform quantum CTRW with different time $`t`$: As shown above in equation (10), quantum CTRW can be described as unitary operator $`U\left(t\right)`$, this is performed with pulse sequence shown in the following (Note that the global phase $`e^{2i\gamma t}`$ of $`U\left(t\right)`$ is safely ignored in our experiments, since $`\rho \left(t\right)=U\left(t\right)\rho \left(0\right)U\left(t\right)^{}`$, this global phase is meaningless and has no effect on the result of experiment)

$$R_x^2\left(\theta \right)R_y^{12}\left(\pi /2\right)\frac{t}{2}R_x^{12}\left(\pi \right)\frac{t}{2}R_y^{12}\left(\pi /2\right).$$

Here $`R_x^2\left(\theta \right)`$ is equal to $`e^{i\frac{\theta }{2}\widehat{\sigma }_x}`$ that act on the second spin,where $`\theta =n\pi /6`$ and $`t=n/\left(6J\right)=n\pi /\left(6\gamma \right)`$ for $`n\{1,2,\mathrm{},12\}`$, $`R_x^{12}\left(\pi \right)`$ denotes $`\pi `$ non-selective pulse that acts on both qubits about $`\widehat{x}`$. It is obviously that the final state $`\rho \left(t\right)`$ of the quantum CTRW prior to measure is given by $`\rho \left(t\right)=U\left(t\right)\rho \left(0\right)U\left(t\right)^{}`$.

Finally, readout the result $`\rho \left(t\right)`$ and calculate quantum total variation distance $`\mathrm{\Delta }^Q\left(t\right)`$: In NMR experiment, it is not practical to determine the final state directly, but an equivalent measurement can be made by so-called quantum state tomography to recover the density matrix $`\rho \left(t\right)=|\psi \left(t\right)\psi \left(t\right)|`$. However, as only the diagonal elements of the final density operators are needed in our experiments, the readout procedure is simplified by applying gradient pulse before readout pulse to cancel the off-diagonal elements. Here we shall note that, the gradient pulse can remove off-diagonal terms since we use heteronuclear systems in our experiment. Then quantum total variation distance $`\mathrm{\Delta }^Q\left(t\right)`$ is determined by the equation $`\mathrm{\Delta }^Q\left(t\right)=\frac{1}{2}_{k=0}^3\left|P_k^Q\left(t\right)\frac{1}{4}\right|`$, where $`P_k^Q\left(t\right)=k\left|\rho \left(t\right)\right|k`$ is certain probability on the node $`|k`$. Finally, $`\mathrm{\Delta }^Q\left(S\right)`$ is determined with equation (12).

All experiments are conducted at room temperature and pressure on Bruker Avance DMX-500 spectrometer in Laboratory of Structure Biology, University of Science and Technology of China. FIG. 3 show the quantum total variation distance $`\mathrm{\Delta }^Q\left(t\right)`$ as a function of time $`t`$ and FIG. 4 show the quantum total variation distance $`\mathrm{\Delta }^Q\left(S\right)`$ as a function of entanglement $`S`$ of $`|\psi \left(t\right)`$ shown in equation (11). From FIG. 3 and FIG. 4, it is clearly seen the good agreement between theory and experiment. However, there exist small errors which increase when time $`t`$ increase, we think that the most errors are primarily due to decoherence, because the time used to implement quantum CTRW $`U\left(t\right)`$ is increased from several to several tens milliseconds approximately, while the decoherence time $`T_20.3`$ and $`0.4`$ seconds for carbon and proton respectively. The other errors are due to inhomogeneity of magnetic field, imperfect pulses, and the variability over time of the measurement process.

## IV Conclusion

We present the experimental implementation of the quantum random walk algorithm on a two-qubit nuclear magnetic resonance quantum computer. For the quantum CTRW on a circle with four nodes, we observe that the quantum walk behaves greatly differently from its classical version. The quantum CTRW can yield an exactly uniform distribution, and is reversible and periodic, while the classical walk is essentially dissipative. Further, we find that the property of this quantum walk strongly depends on the quantum entanglement between the two qubits. The uniform distribution could be obtained only when the two qubits are maximally entangled. In this paper only the relatively simple case with two qubits are considered. However, our scheme could be extended to the case of a graph containing arbitrary $`N`$ nodes, and the quantum random walk could be carried out by using $`\mathrm{log}_2N`$ qubits.

###### Acknowledgements.

We thank T. Rudolph and Y.D. Zhang for useful discussion. Part of the ideas were originated while J.F. Du was visiting Service de Physique Théorique, Université Libre de Bruxelles, Bruxelles. S. Massar and N. Cerf are gratefully acknowledged for their invitation and hospitality. J.F. Du also thanks Dr. J.F. Zhang for the loan of the chloroform sample. This work was supported by the National Nature Science Foundation of China (Grants No. 10075041 and No. 10075044), the National Fundamental Research Program (2001CB309300) and the ASTAR Grant No. 012-104-0040.