text:
# Heat Transport through Rough Channels

## I Introduction

The study of transport phenomena in irregular media is fully recognized nowadays as an important research subject with many technological and industrial applications. In particular, the role of the morphology of an active interface on the efficiency of heat and mass transport devices has also been a theme of great interest for the scientific community. For example, the development of a modern approach for the project of heat exchangers must necessarily take the detailed description of the surface geometry as a possibility for the proper design of these devices. It is well-known in this field that the inclusion of extended surfaces or fins to the surface of the heat exchanger can substantially increase the rate of heat transfer in the system and, therefore, the efficiency of the equipment Arpaci66 ; Incropera90 .

For all practical purposes, the restrictions imposed to the geometrical features of an exchange interface should be the result of an intimate balance between its complexity and cost from one side, and the performance of the equipment from the other. In any case, an important issue in the design of heat exchangers is certainly the extent to which the heat transfer capacity of the system is dictated by the available surface area. For instance, under diffusion-limited conditions, the access of heat to an irregular surface can be highly nonuniform due to screening effects. In this situation, the regions corresponding to fins or extended protrusions will display a higher activity as opposed to the deep parts of fjords which are more difficult to access. As a result, under strong diffusive constraints, the efficiency of the surface for heat exchange can be substantially different from the one expected, considering its total area, the intrinsic conductive properties of the material, and the transport properties of the heating (or cooling) fluid.

An analogous problem can be found in the context of electrochemistry. The nonuniform activity during the linear transport through irregular interfaces of electrodes has been explained in terms of a screening effect which is entirely similar to the one observed in heat transfer. The extensive research developed on this subject in recent years Sapoval94 ; Pfeifer95 ; Sapoval96  has been mainly devoted to the introduction, calculation and application of the concept of active zone in the Laplacian transport to and across irregular interfaces. For example, through the coarse-graining method proposed in Sapoval94  it is possible to compute the flux through electrode surfaces from their geometry alone, without solving the general Laplace problem. In a subsequent study Sapoval01 , it has been shown that this technique provides consistent predictions for the activity of irregular catalyst surfaces. More recently, the activity of irregular absorbing interfaces operating under diffusion-limited conditions has been investigated through nonequilibrium molecular dynamics Andrade03 . The simulation results indicate that the extent of the interface that is significantly active is rather sensitive to the governing mechanism of transport. Precisely, the length of the active zone decreases continuously with density from the Knudsen to the molecular diffusion regime.

The study of laminar-forced-convective heat transfer through tubes or between parallel plates with smooth walls is well known in the literature as the Graetz problem Graetz . Graetz solved the problem analytically with the assumptions of steady, irrotational and incompressible flow, constant fluid properties, fully developed velocity profiles, and in the absence of energy dissipation effects. Much effort has been dedicated in the past years to improve the accuracy of the Graetz solution and generalize the Graetz problem to include other types of boundary conditions Randall97 ; Dennis56 ; Coelho01 ; Telles02 ; Sellars56 ; Brown60 ; Mikhailov97 . The development and solution of more realistic models for heat transport in internal flow configurations and under forced convection conditions has also been a research theme of great scientific interest Johnston94 ; Papoutsakis81 ; Lahjomri03 . Few studies, however, have been dedicated to the investigation of an important extension of the Graetz problem, namely, the case in which the flow and heat transfer processes are confined to walls of irregular geometry.

In the present work we apply the concept of active zone to investigate the steady-state transport of heat in a fluid flowing through a two-dimensional channel whose walls are irregular interfaces. We investigate through direct computational simulations the effect of the interface geometry on the heat exchange efficiency of the system for different diffusive(conductive)-convective conditions. Compared to the behavior of a smooth channel under conditions where the convective mechanism of transport is relevant, our results show that the activity of the interface can be significantly underestimated if geometrical details are not adequately considered.

This paper is organized as follows. In Section 2, we introduce the general description of the physical system and study the flow through the irregular geometry. The transport of heat through the flowing fluid in the rough channel is then numerically treated in Section 3. Finally, general conclusions and perspectives are drawn in Section 4.

## II Viscous flow between rough interfaces

As shown in Fig. 1, our physical system is a two-dimensional channel of length $`L`$ whose delimiting walls are identical interfaces with arbitrary geometry separated by a distance $`2h`$. By definition, the two interfaces are made by the connection of several wall subsets. In the particular roughness model adopted here, we consider that these subsets are identical pre-fractal curves with the geometry of a square Koch tree Evertsz92 .

The mathematical description for the detailed fluid mechanics in the channel is based on the assumptions that we have a continuum, Newtonian and incompressible fluid flowing under steady state conditions. Thus, the Navier-Stokes and continuity equations reduce to

$`\rho \left[u{\displaystyle \frac{u}{x}}+v{\displaystyle \frac{u}{y}}\right]`$ $`=`$ $`{\displaystyle \frac{p}{x}}+\mu \left[{\displaystyle \frac{^2u}{x^2}}+{\displaystyle \frac{^2u}{y^2}}\right],`$ (1)
$`\rho \left[u{\displaystyle \frac{v}{x}}+v{\displaystyle \frac{v}{y}}\right]`$ $`=`$ $`{\displaystyle \frac{p}{y}}+\mu \left[{\displaystyle \frac{^2v}{x^2}}+{\displaystyle \frac{^2v}{y^2}}\right],`$ (2)
$`{\displaystyle \frac{u}{x}}+{\displaystyle \frac{v}{y}}`$ $`=`$ $`0.`$ (3)
Here the independent variables $`x`$ and $`y`$ denote the position in the channel, $`u`$ and $`v`$ are the components of the velocity vector in the $`x`$ and $`y`$ directions, respectively, and $`p`$ is the local pressure in the system. The relevant physical properties of the flowing fluid are the density $`\rho `$ and the viscosity $`\mu `$. In our simulations, we consider non-slip boundary conditions at the entire solid-fluid interface. In addition, the changes in velocity rates are assumed to be zero at the exit $`x=L`$ (gradientless boundary conditions), whereas a parabolic velocity profile is imposed at the inlet of the channel,

$$u(0,y)=\frac{3}{2}V[1(\frac{y}{h})^2].$$
(4)
The Reynolds number is defined here as $`Re\rho Vh/\mu `$, where $`V`$ is the average velocity at the inlet. For simplicity, we restrict our study to the case where the Reynolds number is sufficiently low ($`Re1`$) to ensure a laminar viscous regime for fluid flow. Finally, we also assume that the fluid density and viscosity are constant properties of the flow throughout the channel. As a consequence, the flow field is decoupled from temperature and can be independently calculated.

The numerical solution of Eqs. (1)-(3) for the velocity and pressure fields in the rough channel is obtained through discretization by means of the control volume finite-difference technique Patankar80 . For the complex geometry involved, we solve this problem using a structured mesh based on quadrangular grid elements. For example, in the case of the channel with walls subsets that are 3-generation Koch curves, a total of 79,200 cells generates satisfactory results when compared with numerical meshes of small resolution. The integral form of the governing equations is then considered at each cell of the numerical grid to produce a set of algebraic equations which are pseudo-linearized and solved using the SIMPLE algorithm Patankar80 . The criteria for convergence used in the simulations is defined in terms of residuals, i.e., a measure of the degree to which the conservation equations are satisfied throughout the flow field. In all simulations we performed, convergence is considered to be achieved only when each of the residuals fall below $`10^6`$.

In Fig. 2 we show the velocity profiles corresponding to two different sections of the rough channel along its longitudinal direction. Although very small when compared to the velocities at the center, the velocities of the fluid cells constituting the reentrant zones are still finite. Indeed, as depicted in Fig. 3, a close-up of these roughness details shows fluid layers in the form of consecutive eddies whose intensities fall off in geometric progression Moffat64 ; Leneweit99 . Although much less intense than the mainstream flow, these recirculating structures are located deeper in the system, and therefore experience closer the landscape of the solid-fluid interface. More precisely, viscous momentum is transmitted laterally from the mainstream flow and across successive laminae of fluid to induce vortices inside the fractal cavity. These vortices will then generate other vortices of smaller sizes and intensities. In the next section, we show how this flow structure affects the overall transport of heat as well as the distribution of fluxes among the active wall elements in the rough channel.

## III Heat transport

Once the velocity and pressure fields are obtained for the flow in the rough channel, we proceed with the mathematical modeling and numerical computation of the temperature scalar field. As already mentioned, we assume that the flow is independent of the temperature. In practice, we simply “freeze” the velocity field solution previously calculated and use it to compute the steady-state temperature field $`T(x,y)`$ solving the two-dimensional diffusion-convection equation,

$$u\frac{T}{x}+v\frac{T}{y}=\alpha \left(\frac{^2T}{x^2}+\frac{^2T}{y^2}\right),$$
(5)
where $`\alpha k/\rho c_p`$ is the thermal diffusivity, $`k`$ is the thermal conductivity, and $`c_p`$ is the specific heat of the fluid. For boundary conditions, we consider that the inlet is subjected to a constant temperature $`T_0`$ and that the temperature $`T_w`$ is imposed on both wall interfaces. In particular, we study the case in which $`T_0>T_w`$, i.e., the heat is transported from the bulk of the fluid to the walls of the rough channel. The relative influence of the conductive (diffusive) and convective mechanisms of heat transport on the behavior of the system is quantified here in terms of the Péclet number, defined as

$$PeVh^2/\alpha L.$$
(6)
Similarly to the procedure already utilized for the calculation of the flow, the temperature $`T(x,y)`$ is also obtained through discretization of Eq. (5) by means of finite-differences, with an upwind scheme to avoid numerical instabilities due to the presence of convection Patankar80 . Using this numerical technique, simulations have been performed for three types of channels corresponding to the first three pre-fractal generations of the square Koch tree.

In Fig. 4 we show contour plots of the temperature field in a 3-generation channel for four different diffusion-convection conditions. At very low $`Pe`$ (Fig. 4a), conduction is the dominant mechanism with the transport of heat being practically limited to the entrance region of the channel. By increasing the value of $`Pe`$ (Figs. 4b-d), the contribution of convection to heat transfer through the center of the channel becomes relevant. As shown in Fig. 4c, although the inlet temperature can reach the exit at $`Pe200`$, the temperature at the reentrant zones of the wall remains virtually unaffected by convection. Only at very high values of $`Pe`$, when the vortices inside these irregularities can effectively transport heat, the high temperature front can penetrate deeper in the system (Fig. 4e).

The role of the wall roughness on the transport of heat in the channel can be quantified in terms of the normalized heat flux

$$\varphi q_w/q_0,$$
(7)
where $`q_w`$ is the total heat flux crossing the fluid-solid interface from the bulk,

$$q_w=k_{surf}\frac{T}{n}𝑑s.$$
(8)
The reference value $`q_0`$ corresponds to the heat flux crossing a smooth interface of size $`L`$ at temperature $`T_w`$ that is part of a system where heat is transported by conduction through a static fluid from a line at a fixed temperature $`T_0`$,

$$q_0=k\frac{T_0T_w}{h}L.$$
(9)
The dependence of the flux $`\varphi `$ on the $`Pe`$ number is shown in Fig. 5 for the first three generations of the selected pre-fractal geometry. The behavior of a smooth channel of size $`L`$ and subjected to the same boundary conditions is also shown for comparison. In this simpler case, the system displays a sharp increase of the flux with $`Pe`$, up to a point ($`Pe10`$) where the following power-law regime persists for more than 5 orders of magnitude:

$$\varphi Pe^{1/3}.$$
(10)
This behavior has been observed in previous theoretical studies Sellars56 . In the case of the power-law, it has been demonstrated that the analytical solution of the classical Graetz problem, where the mechanism of axial conduction is disregarded, converges asymptotically to this scaling behavior Sellars56 . Alternatively, as we show in the Appendix, this typical power-law regime of heat transport, known to be valid for moderate and high $`Pe`$ values, can be readily recovered by means of a global heat balance on the system combined with a simple scaling argument.

Remarkably, the results shown in Fig. 5 indicate that different rough channels display approximately the same behavior for low and moderate $`Pe`$ numbers ($`Pe<200`$). Moreover, their performances are also similar to the behavior of the smooth channel for the same range of $`Pe`$ values. Under these conditions, the conduction of heat either predominates over convection or still represents a significant contribution to the overall flux across the walls.

For Péclet values in the range $`10<Pe<200`$, each pre-fractal subset of the interface with the surrounding fluid can be approximately taken as a purely diffusive cell subjected to a constant temperature $`T_w`$ at the wall, and a temperature $`T_0`$ at a line in the bulk of the fluid separated by a given distance from the center of the channel. Considering steady-state conditions, the system can then be treated as a two-dimensional Laplacian transport unit subjected to Dirichlet’s boundary condition. This is a typical situation where the the theorem of Makarov holds Makarov85 . This theorem states that the information dimension of the harmonic measure on a singly connected interface in $`d=2`$ is exactly equal to 1, where $`d`$ represents dimension. Makarov’s theorem has a simple but important consequence: the screening effect due to geometrical irregularities of the interface can be accessed in terms of the ratio $`SL_p/L`$ Sapoval94 , where $`L_p`$ is the perimeter of the interface. Because the active length of the interface $`L_a`$ is of the order of the system size, $`L_aL`$, then $`L_aL_p/S`$ and the factor $`S`$ can be considered as the “screening factor” of the Dirichlet-Laplacian field. When applied to our physical system, the practical meaning of this exact result is that, irrespective of the surface morphology, the length of the zone within the heat exchange wall which receives most of the heat flux should be of the order of the system size $`L`$, instead of the perimeter $`L_p`$. This explains why in Fig. 5 the heat flux at the interface is nearly insensitive to differences in the roughness geometry for low and moderate values of $`Pe`$.

For $`Pe>200`$, the predominance of convection over conduction is extended to the reentrant zones of the pre-fractal units. This crossover value results from the fact that the magnitudes of the velocity vectors constituting the largest vortices in the rough region are approximately two orders of magnitude smaller than the main flow velocity $`V`$. At high $`Pe`$ values, heat can then be convected through these recirculating structures of the flow towards the less accessible wall elements of the interface. Under such circumstances, the theorem of Makarov is not valid anymore, and the flux $`\varphi `$ across the wall becomes dependent on the geometrical details of the interface. As shown in Fig. 5, the higher the generation of the pre-fractal interface, larger is the heat flux at high Péclet conditions.

We now show that the previous analysis based on the overall heat flux $`\varphi `$ is consistent with the behavior of a more detailed measure of transport. In particular, the flux heterogeneity at the interface can be quantified in terms of an active length defined as Sapoval99

$$L_a1/\underset{i=1}{\overset{L_p}{}}\varphi _i^2(1L_aL_p),$$
(11)
where the sum is over the total number of interface elements $`L_p`$, $`\varphi _iq_i/q_j`$, and $`q_i`$ is the thermal flux at the wall element $`i`$. From the definition (11), $`L_a=L_p`$ indicates a limiting state of equal partition of normalized fluxes ($`\varphi _i=1/L_p`$, $`i`$) whereas $`L_a=1`$ should correspond to the case of a fully “localized” flux distribution. The active length therefore provides an useful index to quantify the interplay between the Laplacian phenomenology and the complex geometry of the absorbing interface at the local scale. As expected, the results shown in Fig. 6 indicate that $`L_a`$ generally increases with $`Pe`$. It is important to note, however, that the activities of smooth and rough interfaces behave approximately in the same way for low values of $`Pe`$. This is also consistent with Makarov’s theorem under conditions where the conduction mechanism controls the transport of heat inside the pre-fractal roughness.

For high Péclet values, $`L_a`$ becomes strongly sensitive to the geometrical details of the interface. Again, due to the relevant contribution of convection to the transport at the smaller scales, the flux distribution among the interface units becomes more uniform. The higher the pre-fractal generation used to create the roughness geometry, higher is the available number of wall elements for heat exchange. As a consequence, a substantial increase in the active length can be observed for a fixed value of $`Pe`$.

## IV Conclusions

The development of a coherent scientific approach to the project of heat and mass transport devices involving new concepts from mathematics and physics still represents an important research and technological challenge to chemical and mechanical engineers. In the present work, we have studied the conditions under which the roughness geometry of a channel with a flowing fluid should play a relevant role on the transport of heat. For low Péclet numbers, when the mechanism of heat conduction dominates, it is shown that the overall flux of heat across a pre-fractal interface subjected to Dirichlet boundary conditions has little sensitivity to details and irregularities of the heat absorption zone. Under similar conditions, the length of the active zone computed for a smooth channel is very close to the measure obtained for a rough system. This behavior is then interpreted in terms of a screening phenomenon that privileges the absorbing activity of the geometrically exposed sites of the interface. Moreover, such a “localization” of the heat flux can be formally explained in terms of Makarov’s theorem for Laplacian transport in two-dimensional irregular geometries.

The situation is entirely different at high Péclet numbers. In this case, our simulations reveal that both the heat flux and the active length start to depend noticeably on the morphological details of the interface. In this situation, not only the total heat flux leaving the system increases faster for higher generations of the pre-fractal interface, but the distribution of fluxes at the wall units also becomes gradually less “localized” ($`L_a`$ increases) as $`Pe`$ increases. This transition reflects the onset of convective effects on the heat transport in the proximity of the rough interface. A closer look at the system indicates that the significant changes observed in $`L_a`$ above the transition point are induced by the presence of recirculation zones of flow (vortices) in these regions.

As a future work, we will extend the approach introduced here to simulate the effect of different geometries on the transport efficiency of the system as well as other types of “absorption” mechanisms (e.g., finite-rate heat transfer at the walls) limited by diffusion transport. Finally, we also believe that the results presented in this study provide useful information to the understanding of other physico-chemical systems with relevance to science and technology that include heat, mass or momentum transport. The modeling strategy devised here should also be helpful as a design tool to choose a suitable interface geometry for a given heat exchange system.

## V Acknowledgments

We thank Bernard Sapoval and Marcel Filoche for many helpful discussions. We acknowledge the Brazilian agencies CNPq, CAPES and FUNCAP for financial support.

## VI Appendix: The smooth channel

Here we revisit the classical problem of heat transport through a fluid flowing in a channel whose walls are smooth interfaces. The system setup is shown in Fig. 7. As in the case of the rough channel, we consider a Newtonian and incompressible fluid with constant properties that is flowing under viscous and steady-state conditions. Assuming non-slip boundary conditions at the walls, the velocity along the channel follows the classical parabolic profile

$$u=\frac{3}{2}V(1\eta ^2),$$
(12)
where $`\eta y/h`$, $`V\mathrm{\Delta }ph/(3\mu L)`$, and $`\mathrm{\Delta }p`$ is the pressure drop between the entrance and the exit of the channel. Under these simple flow conditions, the fluid is suddenly subjected to the temperature $`T_0`$ at the channel inlet, while a constant temperature $`T_w`$ is imposed at both walls. At steady-state, and assuming that the axial heat conduction mechanism is negligible, this corresponds to the so-called Graetz problem, where the temperature field in the channel obeys the following differential equation:

$$u(y)\frac{T}{x}=\alpha \frac{^2T}{y^2},$$
(13)
with boundary conditions

$`T(0,y)=T_0\text{and}T(x,h)=T_w.`$ (14)
An analytical solution for the linear partial differential equation (13) subjected to (14) can be obtained through separation of variables and a proper expansion in terms of the so-called Graetz functions Arpaci66 ; Sellars56 . From this solution, it is possible to show that the overall heat flux through the absorbing walls, $`\varphi `$, follows the power law behavior, $`\varphi Pe^{1/3}`$, that is valid for moderate and high $`Pe`$ values.

The objective here is to show that this power-law regime of heat transport can be readily recovered by means of a global heat balance combined with a simple scaling argument. First, in the absence of sources or sinks of heat, the total heat crossing the walls should be equal to the difference

$$q_w=q_{in}q_{out},$$
(15)
where $`q_{in}`$ and $`q_{out}`$ are the fluxes at the entrance and at the exit sections of the channel, respectively. At sufficiently high values of the Péclet number, the convective part of the flux dominates the heat transfer process. As shown in Fig. 8, the temperature profile at the exit displays a central plateau in the range $`y^{}yy^{}`$, where the temperature is equal to the entrance temperature $`T_0`$, and then decays monotonously towards the wall temperature $`T_w`$, in the region $`y^{}|y|h`$. The value of $`\eta ^{}(y^{}/h)`$ can be found from the definition of the local Péclet number

$$Pe_\eta \frac{\tau _D}{\tau _C}=\frac{\delta ^2/\alpha }{L/u(\eta )},$$
(16)
where $`\tau _D\delta ^2/\alpha `$ is the characteristic conduction time, $`\tau _CL/u`$ is the local convective time, and $`\delta h(1\eta ^{})`$. If we now argue that $`\eta ^{}`$ should correspond to the position where $`Pe_\eta 1`$, we can write from Eqs. (12) and (16) that

$$\frac{3}{2}\frac{Vh^2}{\alpha L}(1\eta ^{})^2[1(\eta ^{})^2]=1,$$

which, from the definition (6), leads to

$$(1\eta ^{})^3(1+\eta ^{})=\frac{2}{3Pe}.$$

For high values of $`Pe`$, $`\eta ^{}`$ is very close to 1, such that $`(1\eta ^{})^3(1+\eta ^{})2(1\eta ^{})^3`$, resulting in

$$\eta ^{}1\frac{1}{(3Pe)^{1/3}}.$$
(17)
Going back to Eq. (15), the heat flux $`q_{in}`$ entering the system is given by

$$q_{in}=2\rho c_pT_0_0^hu(y)𝑑y.$$

The flux $`q_{out}`$ leaving the channel can be expressed as the sum

$$q_{out}=q_1+q_2,$$
(18)
where the first term corresponds to the heat escaping through the plateau region of temperature $`T_0`$ and size $`y^{}`$

$$q_1=2\rho c_pT_0_0^y^{}u(y)𝑑y,$$
(19)
while the second term represents the flux through the remaining part of the channel cross-section with variable temperature

$$q_2=2\rho c_p_y^{}^hu(y)T(y)𝑑y.$$
(20)
As a first approximation, if we assume that the temperature in this region decreases linearly with $`y`$,

$$T(y)=\frac{T_0T_w}{y^{}h}(yy^{})+T_0,$$

Eq. (20) can then be written as

$`q_2`$ $`=`$ $`2\rho c_pT_0{\displaystyle _y^{}^h}u(y)𝑑y`$
$`+`$ $`2\rho c_p{\displaystyle _y^{}^h}u(y){\displaystyle \frac{T_0T_w}{y^{}h}}(yy^{})𝑑y,`$
and the total flux becomes

$$q_w=2\rho c_p\frac{T_0T_w}{hy^{}}_y^{}^h\frac{3}{2}V\left[1\left(\frac{y}{h}\right)^2\right](yy^{})𝑑y.$$

Rewriting this expression in terms of the dimensionless variable $`\eta `$, we obtain

$`q_w`$ $`=`$ $`3\rho c_p(T_0T_w)Vh{\displaystyle \frac{1}{1\eta ^{}}}{\displaystyle _\eta ^{}^1}\left(1\eta ^2\right)(\eta \eta ^{})𝑑\eta `$
$`=`$ $`{\displaystyle \frac{3}{12}}\rho c_p(T_0T_w)Vh\left[(\eta ^{})^3+(\eta ^{})^25\eta ^{}+3\right].`$
From this equation and the definitions (6) and (7), the normalized flux can be expressed as

$$\varphi =\frac{3}{12}Pe\left[(\eta ^{})^3+(\eta ^{})^25\eta ^{}+3\right],$$

and using the approximation (17), we finally get the scaling relation

$$\varphi \frac{\left(3Pe\right)^{1/3}}{3}\frac{1}{12}.$$
(21)


format: text

text:
# “Phase Transitions in the Early Universe” and “Defect Formation” 1footnote 11footnote 1Lectures given at a COSLAB Workshop on “Cosmological Phase Transitions and Topological Defects”, 22–24 May 2003, Porto, Portugal.

## I Introduction

If one calculates the equilibrium properties of the Standard Model at high temperatures, one finds several different phases Kapusta ; Linde ; Bailin . Below 100 MeV, the elementary particles behave more or less as at zero temperature. If we heat them up above this temperature, there is a transition to a deconfined phase in which quarks and gluons are not confined into hadrons. At 100 GeV, there is another transition at which the weak and electromagnetic interactions unify. Extensions of the Standard Model, such as Grand Unified Theories (GUTs), typically predict new phases at even higher temperatures.

This phase structure is important in cosmology Kolb ; Linde  because the temperature of the universe was initially very high. It seems likely that had to be higher than 100 GeV, because otherwise it is extremely difficult to explain the baryon asymmetry of the universe — the fact that there is so much matter in the universe, but practically no antimatter. As the universe expanded and cooled down to the current 2.7 K, it must have therefore undergone at least the deconfinement and electroweak phase transitions mentioned above. Typical inflationary models predict much higher temperatures, usually at least $`10^{10}`$ GeV, and it is therefore likely that there were also other phase transitions.

The cosmological consequences of these phase transitions depend sensitively on their properties. In many cases the consequences would be disastrous, and one can then rule out the theory that predicts the transition, or at least constrain its parameters. Examples of this are grand unified theories with strongly first-order phase transitions, as they would lead to a period of “old inflation”, and to a very inhomogeneous universe, or theories that predict massive domain walls or too many magnetic monopoles. In other cases, the effects of the phase transition may be subtle, only barely observable, such as gravitational waves, or even completely unobservable. Or, it may be that the consequences are absolutely crucial for our own existence, as in the theory of electroweak baryogenesis Shaposhnikov , which uses the electroweak phase transition to explain the baryon asymmetry.

In this lecture, I will concentrate on those phase transitions that can be understood as a spontaneous breakdown of either a global or a local symmetry. This covers the electroweak and GUT transitions, and to a certain extent also the deconfinement transition. I will also point out similarities to phase transitions in condensed matter systems such as superfluids and superconductors. Unlike cosmological phase transitions, they can be studied experimentally and can therefore be used to test our theoretical understanding of phase transitions in quantum field theories.

## II Definitions

The properties of a thermodynamic system in thermal equilibrium are given by the partition function $`Z`$ (see Ref. Kapusta ), which can be generally written as

$$Z=\mathrm{Tr}\widehat{\rho }=\mathrm{Tr}e^{\beta \widehat{H}},$$
(1)
where $`\widehat{H}`$ is the Hamiltonian of the system and $`\beta =1/kT`$ is the inverse temperature.

In many cases, it is more convenient to use its logarithm, which defines the free energy $`F=T\mathrm{ln}Z`$. This quantity has the property that its value is minimised in thermal equilibrium.

Neither the partition function nor the free energy are observable quantities. However, any observable can be expressed as the derivative of $`F`$ with respect to some external field. For instance, the expectation value of the operator $`\widehat{X}`$ is defined as

$$\widehat{X}=Z^1\mathrm{Tr}\widehat{X}\widehat{\rho },$$
(2)
and if we define an external field $`J_X`$ by

$$Z(J_X)=\mathrm{Tr}e^{\beta \widehat{H}+J_X\widehat{X}},$$
(3)
we find

$$\widehat{X}=\frac{}{J_X}\mathrm{ln}Z=\frac{1}{T}\frac{F}{J_X}.$$
(4)
A phase transition is defined as a point at which the free energy (or, equivalently, the partition function) is non-analytic. In general, it follows that at least some observables are non-analytic at the transition point and act as order parameters.

The simplest type of phase transitions are first-order transitions, in which some expectation value, say $`\widehat{X}`$, is discontinuous. An everyday example of this is the transition between water and ice. Let us assume that the transition takes place at temperature $`T_c`$. Because of the discontinuity, the state of the system at $`T_c`$ depends on which side the critical temperature was approached from. Consequently, it is possible to have two phases coexisting at $`T_c`$, which is exactly what happens when we put ice into a glass of water. As the temperature is gradually changed through $`T_c`$, the discontinuous jump from one value of $`\widehat{X}`$ to the other does not necessarily take place instantaneously at $`T_c`$. For instance, water can remain liquid in a supercooled state at temperatures below $`0^{}`$C for a relatively long time, if it is pure enough and completely unperturbed. This happens in a meteorological condition known as freezing rain, in which the supercooled water freezes as soon as it hits the ground, causing hazardous driving conditions. More generally, this phenomenon is known as metastability.

Second-order phase transitions are those in which every expectation value is continuous, but some expectation value $`\widehat{X}`$ has a discontinuous derivative. This is a limiting case of a first-order transition in which the discontinuity is taken to zero, and therefore it typically requires fine tuning of some parameter. On a phase diagram, second order transitions are typically end points of first order transition lines. One example of this is the critical point of water at which the first-order transition line between water and vapour ends. To reach this point, one has to tune both the temperature and the pressure to $`T=647.096`$ K and $`p=22.064`$ MPa. In contrast, reaching a first order transition point, one only has to tune the temperature provided that the pressure is below the critical value.

There is, however, an exception to this rule. There are cases in which some symmetry forces the system to be on the first-order transition line. This happens, for instance, in a ferromagnet. If the temperature is below the Curie temperature, the magnetisation $`\stackrel{}{M}`$ becomes non-zero even if there is no external magnetic field $`\stackrel{}{H}`$. Because the setting is rotation invariant, no direction of $`\stackrel{}{M}`$ is preferred. Therefore, the direction of $`\stackrel{}{M}`$, which depends on how the external field was taken to zero, breaks the rotation invariance. This is known as spontaneous symmetry breaking.

If we imagine changing the external field $`\stackrel{}{H}`$ gradually through zero, we see that $`\stackrel{}{M}`$ is discontinuous at $`\stackrel{}{H}=0`$. In other words, there is a first-order transition at $`\stackrel{}{H}=0`$. If we increase the temperature keeping $`\stackrel{}{H}`$ zero, we are therefore guaranteed to reach the second-order end point at the Curie temperature.

The above situation is very important in field theories, because they typically have a high degree of symmetry. As in the case of the ferromagnet, it often happens at low temperatures that the state of the system is not invariant under a symmetry of the theory and this typically leads to a second-order phase transition as the temperature is increased.

## III Global symmetry breaking in field theories

As a simple example of a field theory with spontaneous symmetry breaking, let us consider the Lagrangian

$$=\frac{1}{2}_\mu \varphi ^\mu \varphi \frac{\lambda }{4}\left(\varphi ^2v^2\right)^2,$$
(5)
where $`\varphi `$ is a real scalar field. The Lagrangian is symmetric under $`\varphi \varphi `$, which correspond to the group $`_2`$.

The vacuum of the system is found by solving the Euler-Lagrange equations for a constant field,

$$\frac{}{\varphi }=0,$$
(6)
which gives

$$\lambda \left(\varphi ^3v^2\varphi \right)=0.$$
(7)
In addition to $`\varphi =0`$, which is an unstable solution, we find two solutions, $`\varphi =\pm v`$, which are the possible vacuum states of the system. Neither choice is invariant under the $`_2`$ symmetry, and therefore the symmetry is spontaneously broken in the vacuum.

At a non-zero temperature, the value of $`\varphi `$ fluctuates from point to point. If these thermal fluctuations are small, we can write the field as the sum of the homogeneous background contribution $`\varphi _0`$ and the fluctuation $`\delta \varphi `$,

$$\varphi (\stackrel{}{x})=\varphi _0+\delta \varphi (\stackrel{}{x}).$$
(8)
We assume that these fluctuations are Gaussian, which is equivalent to the one-loop approximation in perturbation theory, i.e., in the expansion in terms of Feynman diagrams. In many cases, this is a good approximation. The homogeneous mode $`\varphi _0`$ is chosen in such a way that the mean value of the fluctuation in zero, $`\delta \varphi =0`$. The variance, however, is non-zero, and has the value Kapusta

$$\delta \varphi ^2=\frac{d^3p}{(2\pi )^3}\frac{1}{\omega }\frac{1}{e^{\beta \omega }1}+(\text{zero-temperature part}),$$
(9)
where $`\omega =(\stackrel{}{p}^2+m^2)^{1/2}`$. Here $`m`$ is the mass of the field in the vacuum around which the expansion is done, and is given by the second derivative of the potential,

$$m^2=\frac{^2V}{\varphi ^2}.$$
(10)
The zero-temperature part is actually divergent, but we do not have to worry about that, because it is already absorbed in the definition of the renormalised zero-temperature couplings. For $`Tm`$, Eq. (9) can be expanded as

$$\delta \varphi ^2=\frac{T^2}{12}\frac{mT}{4\pi }.$$
(11)
When the fluctuations $`\delta \varphi `$ are taken into account, the homogeneous more $`\varphi _0`$ feels an effective potential that is different from the “tree-level” potential in Eq. (5). This is easiest to see by averaging Eq. (7) over fluctuations Linde ,

$`0`$ $`=`$ $`\lambda \left(\varphi ^3v^2\varphi \right)`$ (12)
$`=`$ $`\lambda \left(\varphi _0^3+3\delta \varphi ^2\varphi _0v^2\varphi _0\right),`$
where we have used the fact that the fluctuations are symmetric and odd moments must therefore vanish. If we now use the leading term from Eq. (11), $`\delta \varphi ^2T^2/12`$, we find

$$\lambda \left(\varphi _0^2+\frac{T^2}{4}v^2\right)\varphi _0=0.$$
(13)
This shows that the two non-zero minima only exist below a critical temperature $`T_c=2v`$, and above that, the symmetry is restored.

One can easily generalise this discussion to somewhat more complicated cases. If $`\varphi _i`$ is a multicomponent real scalar field with $`i1,\mathrm{},N`$, we can consider the Lagrangian

$$=\frac{1}{2}_\mu \varphi _i^\mu \varphi _i\frac{\lambda }{4}\left(\varphi _i\varphi _iv^2\right)^2,$$
(14)
where summation over the indices $`i`$ is assumed. This Lagrangian is invariant under SO($`N`$) rotations of the field $`\varphi _i`$. For two components, we can define a complex scalar field $`\mathrm{\Phi }=(\varphi _1+i\varphi _2)`$, and we have

$$=\frac{1}{2}_\mu \varphi ^{}^\mu \varphi \frac{\lambda }{4}\left(\varphi ^{}\varphi v^2\right)^2.$$
(15)
This Lagrangian is invariant under rotations of the complex phase angle of the scalar field $`\varphi `$., which means that it has a U(1) symmetry. Because of this symmetry, we can choose that the vacuum expectation value is real, and expand $`\varphi _1=\varphi _0+\delta \varphi _1`$ and $`\varphi _2=\delta \varphi _2`$. The variance of both $`\delta \varphi _1`$ and $`\delta \varphi _2`$ are given by Eq. (11). The averaged Euler-Lagrange equation for the real part is

$`0`$ $`=`$ $`{\displaystyle \frac{}{\varphi _1}}`$ (16)
$`=`$ $`\lambda \left(\varphi _1^3+\varphi _2^2\varphi _1v^2\varphi _1\right)`$
$`=`$ $`\lambda \left(\varphi _0^2+3\delta \varphi _1^2+\delta \varphi _2^2v^2\right)\varphi _0.`$
Using Eq. (11), we find

$$\lambda \left(\varphi _0^2+\frac{T^2}{3}v^2\right)\varphi _0=0,$$
(17)
showing that this time the critical temperature is $`T_c=\sqrt{3}v`$, but otherwise the conclusions are unchanged.

The restoration of broken symmetries at high temperatures is important for the early universe, because the universe was initially extremely hot and because the theories of particle physics contain several symmetries that are spontaneously broken at zero temperature. As shown in Table 1, these are associated with the hypothetical Grand Unified Theory (GUT), the electroweak theory and quantum chromodynamics (QCD).

The QCD deconfinement phase transition takes place at around 100 MeV, and is the only one of the three transitions that can be studied experimentally. The transition cannot be fully understood within perturbation theory, and it can be characterised in two different ways, either as a transition from a confining to a deconfined phase, or as spontaneous breakdown of the chiral symmetry. We will adopt the latter point of view.

If the up and down quarks were massless, the left and right handed quarks would decouple. One could therefore carry out independent rotations between the up and down quarks of the given handedness. In the QCD vacuum, this symmetry is spontaneously broken by the quark condensate $`\overline{q}q`$, which shows how the left and right handed quarks are aligned in the flavour space. Just like in a ferromagnet, one would expect the symmetry to be restored at high temperatures, and this seems to imply a second-order phase transition.

However, the quarks are not actually massless and therefore the chiral symmetry is not exact. The quark mass term plays the same role as the magnetic field $`\stackrel{}{H}`$ in the ferromagnet phase diagram in Fig. 2, and shifts the system a little bit away from the first-order transition line. Because there is no symmetry in the first place, there is no second order transition either. A first-order transition is still possible if the phase diagram is more complicated than that in Fig. 2.

The GUT and electroweak phase transitions, on the other hand, involve exact symmetries, but local gauge symmetries instead of global ones. Therefore, the above discussion does not apply them either, as we shall see in the next section.

## IV Gauge symmetry breaking

The Lagrangians discussed above had global symmetries. One has to carry out the same rotation at every point in order for the Lagrangian to remain unchanged. It is possible to modify the theory in such a way that it is invariant under position-dependent rotations as well. This is done by introducing a gauge field $`A_\mu `$, and the local symmetries are therefore also known as gauge symmetries. For the theory with a complex scalar in Eq. 15), the modified Lagrangian is

$$=\frac{1}{4}F_{\mu \nu }F^{\mu \nu }+\frac{1}{2}D_\mu \varphi ^{}D^\mu \varphi \frac{\lambda }{4}\left(\varphi ^{}\varphi v^2\right)^2,$$
(18)
where we have used the covariant derivative $`D_\mu =_\mu +igA_\mu `$ and the field strength tensor $`F_{\mu \nu }=_\mu A_\nu _\nu A_\mu `$. This model is essentially a relativistic version of the Ginzburg-Landau theory of superconductivity. Our results will therefore apply to superconductors on a qualitative level. The quantitative results differ, though.

The gauge coupling constant $`g`$ is simply the electric charge of the field $`\varphi `$. (For superconductors, $`g=2e`$, because a Cooper pair consists of two electrons.) The Lagrangian is invariant under local (position-dependent) rotations of the phase angle

$$\varphi e^{i\alpha }\varphi ,A_\mu A_\mu \frac{1}{g}_\mu \alpha .$$
(19)
Again, we write $`\varphi =\varphi _0+\delta \varphi `$, where we now choose $`\varphi _0`$ to be real and positive. The gauge field $`A_\mu `$ does not have a homogeneous mode, and therefore we can treat is as a pure fluctuation. Again, we average the Euler-Lagrange equation for the real part $`\varphi _1`$ over fluctuations $`\delta \varphi _1`$, $`\delta \varphi _2`$ and $`A_\mu `$,

$`0`$ $`=`$ $`{\displaystyle \frac{}{\varphi _1}}`$ (20)
$`=`$ $`\lambda \left(\varphi _1^3+\varphi _2^2\varphi _1v^2\varphi _1\right)+g^2A_\mu A^\mu \varphi _1`$
$`=`$ $`\lambda \left[\varphi _0^2+3\delta \varphi _1^2+\delta \varphi _2^2+{\displaystyle \frac{g^2}{\lambda }}A_\mu A^\mu v^2\right]\varphi _0.`$
The gauge field contains three physical degrees of freedom, and therefore we have

$$A_\mu A^\mu =3\delta \varphi ^2=\frac{T^2}{4}\frac{3m_\gamma T}{4\pi }+(T=0\text{ part}).$$
(21)
Crucially, however, the mass $`m_\gamma `$ of the gauge field depends on $`\varphi _0`$ because of the Higgs mechanism, and is $`m_\gamma =g\varphi _0`$. Including this contribution, we find

$$\lambda \left[\varphi _0^2\frac{3}{4\pi }\frac{g^3}{\lambda }T\varphi _0+\left(\frac{1}{3}+\frac{g^2}{4\lambda }\right)T^2v^2\right]\varphi _0=0.$$
(22)
To interpret this results, we define the effective potential as the tree-level potential that would give this same Euler-Lagrange equation. That is obtained by simply integrating the left hand side of the equation,

$`V_{\mathrm{eff}}(\varphi _0)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\left({\displaystyle \frac{\lambda }{3}}+{\displaystyle \frac{g^2}{4}}\right)T^2\lambda v^2\right]\varphi _0^2`$ (23)
$`{\displaystyle \frac{g^3T}{4\pi }}\varphi _0^3+{\displaystyle \frac{\lambda }{4}}\varphi _0^4.`$
The presence of the cubic term means that as long as the coefficient of the quadratic term is positive and small enough, there are two non-negative minima, one at $`\varphi _0=0`$ and the other at a non-zero value. The two minima are degenerate at a certain critical value of $`T_c`$,

$$T_c=\sqrt{\frac{12\lambda v^2}{4\lambda +3g^2(3/2\pi ^2)(g^6/\lambda )}}.$$
(24)
This corresponds to a first-order transition line. Thus we see that due to thermal fluctuations, we have a first-order transition instead of a second-order one.222Note that this is the reason why Type I superconductors have a first-order transition. The discontinuity in $`\varphi _0`$ is

$$\mathrm{\Delta }\varphi _0=\frac{1}{2\pi }\frac{g^3}{\lambda }T_c.$$
(25)
It is also straightforward to calculate the range of temperatures in which a second, metastable minimum exists,

$$\frac{12\lambda v^2}{4\lambda +3g^2}<T^2<\frac{12\lambda v^2}{4\lambda +3g^2(27/16\pi ^2)(g^6/\lambda )}$$
(26)
The critical temperature $`T_c`$, together with the metastability range is shown in Fig. 3 for $`g=0.5`$. Fig. 4 shows the shape of the effective potential at $`T_c`$ and at the ends of the metastability range for $`g=0.5`$, $`\lambda =0.05`$.

When deriving Eq. (23), we have assumed that the fluctuations are Gaussian. It is important to check whether this approximation, which is equivalent to a one-loop calculation in terms of Feynman diagrams, is consistent. At higher orders in the loop expansion, one finds that each loop gives typically a contribution of the order $`g^2T/m_\gamma `$. At $`T_c`$, the photon mass in the broken phase is

$$m_\gamma =g\mathrm{\Delta }\varphi _0=\frac{1}{2\pi }\frac{g^4}{\lambda }T_c,$$
(27)
and therefore each loop order gives a contribution proportional to $`\lambda /g^2`$. If $`\lambda g^2`$, this is small and means that the contribution from higher loop diagrams can be ignored. One can therefore trust the calculation of the effective potential in Eq. (23). On the other hand, if $`\lambda g^2`$, the higher loop orders are not suppressed relative to the one-loop contribution, and the whole approximation breaks down.

In the electroweak theory the ratio $`\lambda /g^2`$ is essentially given by the ratio of the Higgs and $`W`$ boson masses. The Higgs boson has not been found and its mass is therefore unknown, but there is an experimental limit that its mass has to be at least 115 GeV. This means that $`\lambda >g^2`$, and the perturbative calculation breaks down. This is also the case in Type II superconductors. In both cases, one needs non-perturbative methods to study the properties of the phase transition.

## V Non-perturbative results

In the global theory, one can be certain that no matter what non-perturbative effects there may be near the transition, the transition must exist, because one can measure the value of $`\varphi =\varphi _0`$. It vanishes in the symmetric phase and is non-zero in the broken phase, and a function that does that cannot be analytic. In other words, $`\varphi `$ acts as an order parameter.

When one considers non-perturbative effects in the gauge theory, one makes the simple but rather striking observation that $`\varphi `$ cannot be an order parameter, because it is not physically observable Elitzur . This is because it transforms non-trivially under gauge transformations,

$$\varphi (x)e^{i\alpha (x)}\varphi (x).$$
(28)
If one averages over all gauge-equivalent configurations, the expectation value of any such quantity vanishes.

The most reliable way of studying phase transitions non-perturbatively is to use lattice Monte Carlo simulations Smit , and they show that the theory in Eq. (18) still has a transition. At small $`\lambda /g^2`$, the transition is of first order as our perturbative calculation indicated, but when $`\lambda /g^2`$ is larger, it becomes a continuous second-order transition. Even though one cannot use $`\varphi `$ as an order parameter, one can see which phase the system is in by studying the properties of the magnetic field. It is massless in the symmetric phase and massive in the broken phase Peisa . The latter is equivalent to the Meissner effect in superconductors.

In non-Abelian theories, there are generally no massless fields, because the theories are confining. Therefore there is no qualitative difference between the “symmetric” and “broken” phases of the electroweak theory. Careful lattice Monte Carlo simulations have shown that when the Higgs mass is above 80 GeV, the two phases are actually smoothly connected and there is no phase transition between them Kajantie . Because the true Higgs mass has to be in this range, it seems likely that there was never a true electroweak phase transition in the early universe. On the other hand, this conclusion may still change if the Higgs sector of the electroweak theory is more complicated than in the simplest version of the theory.

Similarly, the GUT phase transition may have been a smooth crossover from one phase to another, provided that $`\lambda /g^2`$ was large. This is perhaps surprising, because the broken phase of any GUT is supposed to be the SU(3)$`\times `$SU(2)$`\times `$U(1) phase of the Standard Model, and it contains a massless photon. However, the photon mass is not exactly zero in GUTs Polyakov . It is finite, but suppressed by $`\mathrm{exp}(T_{\mathrm{GUT}}/T)`$, and therefore the photon is massless for all practical purposes. Nevertheless, an arbitrarily small mass is enough to make a smooth crossover possible.

## VI Conclusions

We have seen that symmetries that are broken in the vacuum are typically restored at high temperatures, and that this leads to phase transitions. In global theories, the transitions are generally of second order, and in gauge theories they are of first order if the scalar self coupling $`\lambda `$ is small relative to the gauge coupling $`g^2`$. However, in the opposite limit, the transition may become a smooth crossover between two phases, analogously to the transition between water and vapour at very high pressures.

At first sight, it seems there have been at least three phase transitions related to particle physics in the early universe: the deconfinement transition of QCD, the electroweak transition and the GUT transition. On closer scrutiny it turned out that it is likely that they may have all been merely smooth crossovers rather than true phase transitions in the mathematically rigorous sense. Nevertheless, they may still have important consequences especially in the form of defect formation.

Lecture II:
Defect Formation

## VII Introduction

Topological defects are non-linear, time-independent solutions of field theories, which are made stable by the topology of the theory. They have been studied extensively in cosmology Vilenkin , condensed matter physics Volovik  and, at a more theoretical level, in the context of quantum or classical field theory.

In cosmology, magnetic monopoles formed at the GUT phase transition would have disastrous effects later on in the evolution of the universe, and part of the motivation behind the inflationary theory Guth:1980zm  was to wipe them out before that. Cosmic strings, on the other hand, are fully compatible with what is known about cosmology. If they were formed at some early stage, we may well be able to detect them with astronomical observations.

Topological defects also exist in many different condensed matter systems. The ones that are closest to defects of cosmological importance are vortices in superfluids and Abrikosov flux lines in superconductors.

It was pointed out by Kibble in 1976 Kibble  that topological defects are generally formed at symmetry breaking phase transitions and that we should therefore expect that they were formed in the early universe, too. However, as discussed in Ref. proc1 , the known cosmological phase transitions are not simple symmetry breaking transitions, but they involve a breakdown of a local gauge symmetry. Nevertheless, as we shall see, the same conclusion still holds: Topological defects are formed.

In recent years, defect formation has been studied extensively in condensed matter systems such as liquid crystals TurokNature ; Bowick:1994rz ; Digal:1998ak , superfluids Hendry ; Dodd ; Ruutu:1996qz ; Bauerle  and superconductors Carmi ; Monaco ; Kirtley ; Maniv . The first two of these involve global symmetries, and indeed, they seem to confirm Kibble’s theory. Superconductors have a U(1) gauge symmetry and would therefore be closer to the cosmological phase transitions. The superconductor experiments carried out so far have not produced conclusive results, but hopefully this will change soon.

In this lecture, I will first explain how topological defects form in transitions associated with a breakdown of a global symmetry. This is known as the Kibble-Zurek mechanism. Then, I will discuss gauge symmetries and show that the mechanism is actually rather different in that case. For a more comprehensive review, see Refs. IJMPA ; contphys .

## VIII Discrete global symmetries

Consider first a real scalar field in a potential with a $`_2`$ symmetry,

$$V(\varphi )=\frac{1}{2}m^2\varphi ^2+\frac{1}{4}\lambda \varphi ^4.$$
(29)
If $`m^2`$ is positive, the only minimum is at $`\varphi =0`$, and the state of the system is symmetries under $`_2`$. However, for $`m^2<0`$, there are two degenerate minima at $`\varphi =\pm \sqrt{m^2/\lambda }`$.

If the field is initially exactly zero, it remains zero even if $`m^2`$ is negative. However, in real systems there are always both quantum and thermal fluctuations. Although the arguments we will use in this section are probably also valid for pure quantum fluctuations, we assume that the temperature is non-zero and the dominant fluctuations are thermal. Therefore, we will not consider any quantum effects.

If either the value of $`\varphi `$ or its velocity is initially non-zero, it determines which minimum the system ends up in. If the initial condition is inhomogeneous, i.e., position-dependent, different parts of the system may end up in different vacua. On the whole, each minimum is equally probable, and therefore the system ends up in a configuration of positive and negative regions, which both cover half of the whole system. The boundaries between these domains are topological defects known as domain walls. The system is in vacuum on both sides of the wall, but because of continuity, the field must leave the vacuum inside the wall to interpolate between the positive and negative values. This means that the wall has positive energy per unit area.

The most obvious question about the domain configuration is whether it has a characteristic length scale. In general, the initial fluctuations have a finite correlation length $`\xi `$, which is defined by

$$\varphi (\stackrel{}{x})\varphi (\stackrel{}{y})\varphi (\stackrel{}{x})\varphi (\stackrel{}{y})e^{|\stackrel{}{x}\stackrel{}{y}|/\xi }.$$
(30)
This means that if the points $`\stackrel{}{x}`$ and $`\stackrel{}{y}`$ are further away from each other than $`\xi `$, the initial conditions at these points are independent. If we then assume that the domains form in a short enough time that no signal can travel from $`\stackrel{}{x}`$ to $`\stackrel{}{y}`$ before they have formed, there is a 50% chance that the system will end up in different vacua at these two points, in which case there would have to be (at least) one domain wall between them. Thus, one expects the typical distance between the walls to be roughly $`\xi `$.

Let us now assume that $`m^2`$ is initially positive but decreases with time. Initially, the system is in thermal equilibrium in the symmetric phase, but when $`m^2`$ reaches a critical value, a symmetry breaking phase transition takes place. Because the leading effect of a non-zero temperature is to change the effective mass term proc1 , we can equally well assume that $`m^2`$ is constant, but $`T`$ decreases due to, say, the expansion of the universe or a contact with a colder heat bath.

It is a general property of second-order phase transitions that the correlation length diverges at the transition point. This divergence can be characterised by a critical exponent $`\nu `$, which is universal in the sense that it only depends on the universality class and not on any microscopic details,

$$\xi (T)=\xi _0\left(\frac{T_cT}{T_c}\right)^\nu .$$
(31)
Typical values of $`\nu `$ are around 2/3.

If we could cool the system adiabatically so that it stays in equilibrium, the correlation length $`\xi `$ would be infinite at the transition temperature $`T_c`$. Consequently, the whole system could end up in the same vacuum and there would be no domain walls.

However, this is not possible in practice, because of critical slowing down Zurek:qw . The relaxation time $`\tau `$, i.e., the time it takes for a small perturbation to equilibrate, also diverges exponentially, with exponent $`\mu `$,

$$\tau (T)=\tau _0\left(\frac{T_cT}{T_c}\right)^\mu .$$
(32)
This means that if we are cooling the system with a constant rate characterised by the “quench time” $`\tau _Q`$,

$$T(t)=\left(1\frac{t}{\tau _Q}\right)T_c,$$
(33)
the equilibrium correlation length and relaxation time diverge as

$$\xi (t)=\xi _0\left(\frac{|t|}{\tau _Q}\right)^\nu ,\tau (t)=\tau _0\left(\frac{|t|}{\tau _Q}\right)^\mu .$$
(34)
Once $`\tau (t)>|t|`$, any deviation from equilibrium will not have time to disappear before the transition. The fluctuations in the system at that time freeze out and survive all the way to the transition. In particular this means that the correlation length cannot grow significantly after that time, because it is simply a property of these fluctuations.

This means that the actual correlation length at the time of the transition is

$$\widehat{\xi }=\xi (\widehat{t}),$$
(35)
where $`\widehat{t}`$ is the freeze-out time determined by the equation $`\tau (\widehat{t})=|\widehat{t}|`$. This is straightforward to solve and gives

$$\widehat{t}=(\tau _0\tau _Q^\mu )^{1/(1+\mu )},$$
(36)
whereby the frozen-out correlation length is

$$\widehat{\xi }=\xi _0\left(\frac{\tau _Q}{\tau _0}\right)^{\nu /(1+\mu )}\tau _Q^{\nu /(1+\mu )}.$$
(37)
As we saw before, this is the typical distance between domain walls formed in the transition.

## IX Continuous global symmetries

It is straightforward to generalise these arguments to systems with continuous symmetries. In the case of a complex scalar field with a U(1) symmetry, the set of possible vacua, i.e., the vacuum manifold, is a circle. The vacua are characterised by the phase angle $`\phi `$ of the field $`\varphi =ve^{i\phi }`$. After the transition, the angle $`\phi `$ is uncorrelated at distances longer than $`\widehat{\xi }`$.

Consider now three points A, B and C separated by more than $`\widehat{\xi }`$ (see Fig. 5). Each one has a random value of $`\phi `$, but field $`\varphi `$ is still continuous. It will therefore have to somehow interpolate between the three points, and if the points are not separated by much more than $`\widehat{\xi }`$, it is most likely that the field follows the shortest path on the vacuum manifold, because that costs least energy. This assumption is known as the geodesic rule. However, there is a finite probability that if the geodesic rule is satisfied between A and B, between B and C and between C and A, the phase angle changes by $`2\pi `$ as one follows it from A to B to C to A. In that case, continuity tell us that the field has to leave the vacuum manifold and vanish at some point within the circle formed by these points. That point is a vortex.

As with domain walls, the typical distance between vortices after the transition is $`\widehat{\xi }`$, and their number density is therefore

$$n\widehat{\xi }^2\tau _Q^{2\nu /(1+\mu )}.$$
(38)
It is straightforward to generalise these arguments to more complicated symmetry groups. In particular, one finds that the number density of magnetic monopoles is $`n\widehat{\xi }^3`$.

In principle, Eq. (38) makes a prediction that can be tested in condensed matter experiments. However, the arguments that led to it neglected factors or order one in many places, and therefore even if one knew exactly the values of the parameters $`\xi _0`$ and $`\tau _0`$, the prediction could still easily be off by a factor of ten. A more robust prediction is the dependence on $`\tau _Q`$, but that would require carrying out experiments with a wide range of different cooling rates, and in many cases it is not possible to vary the cooling rate at all.

An even more robust test of the Kibble-Zurek mechanism is given by the spatial distribution of vortices. Let us define $`N_W(R)`$ as the winding number around a circle of radius $`R`$,

$$N_W(R)=\frac{1}{2\pi }_{|\stackrel{}{r}|=R}𝑑\stackrel{}{r}\stackrel{}{}\phi .$$
(39)
On average, this of course vanishes, but its variance $`N_W(R)^2`$ is typically non-zero. After a phase transition in a global theory, the circle passes through around $`2\pi R/\widehat{\xi }`$ domains that are uncorrelated with each other. Every time, the curve moves from one domain to another, the phase angle changes by some random amount $`\mathrm{\Delta }\phi `$. Labelling the domains by $`i`$, assuming the geodesic rule, and ignoring all factors of $`2\pi `$, the total winding number is then

$$N_W(R)\underset{i=1}{\overset{R/\widehat{\xi }}{}}\mathrm{\Delta }\phi _i.$$
(40)
Because the changes $`\mathrm{\Delta }\phi _i`$ are uncorrelated, the total variance is

$`N_W(R)^2`$ $``$ $`{\displaystyle \underset{i,j}{}}\mathrm{\Delta }\phi _i\mathrm{\Delta }\phi _j={\displaystyle \underset{i}{}}\mathrm{\Delta }\phi _i^2`$ (41)
$`=`$ $`\mathrm{\Delta }\phi ^2R/\widehat{\xi }R/\widehat{\xi }.`$
More generally, one can define a scaling exponent $`\nu _D`$ by Digal

$$N_W(R)^2R^{4\nu _D}.$$
(42)
The above result then implies that the Kibble-Zurek mechanism predicts $`\nu _D=1/4`$. This can be compared with a fully random distribution, in which the sign $`s_a`$ of each vortex $`a\{1,\mathrm{},N_{\mathrm{vort}}\}`$ inside the curve is independent of all the others. In that case, the winding number is

$$N_W(R)=\underset{a=1}{\overset{N_{\mathrm{vort}}}{}}s_a,$$
(43)
and consequently

$`N_W(R)^2`$ $``$ $`{\displaystyle \underset{a,b}{}}s_as_b={\displaystyle \underset{a}{}}s_a^2`$ (44)
$`=`$ $`N_{\mathrm{vort}}R^2.`$
Thus, a random distribution of vortices would have $`\nu _D=1/2`$.

## X Gauge symmetries

In a gauge field theory, one cannot apply the Kibble-Zurek arguments as such, because the phase angle $`\phi `$ is not a gauge invariant quantity. Furthermore, while a global system tries to minimise phase gradients $`\stackrel{}{}\phi `$ in order to minimise energy, the condition for minimum energy in a gauge theory is $`\stackrel{}{}\phi +e\stackrel{}{A}=0`$. If the gauge field $`\stackrel{}{A}`$ is non zero, the phase $`\phi `$ does not even attempt to become uniform.

In the special case in which the magnetic field is zero everywhere, one can choose a gauge in which $`\stackrel{}{A}=0`$. Then, the gradients of $`\phi `$ cost energy, and one can expect the same arguments to hold as in a global theory.

However, if one considers a non-zero temperature, as we have been doing, one cannot make this assumption, because it is not consistent with thermal fluctuations Hindmarsh:2000kd . In fact, these fluctuations are nothing but thermal radiation, and in the symmetric phase, they can be approximated by the blackbody spectrum. To the extent that we can consider classical dynamics, we can therefore say that the initial conditions for the magnetic field are given by the Rayleigh-Jeans spectrum. A convenient way to express that is by defining the function $`G(k)`$,

$$B_i(\stackrel{}{k})B_j(\stackrel{}{k}^{})=G(k)\left(\delta _{ij}\frac{k_ik_j}{k^2}\right)(2\pi )^3\delta (\stackrel{}{k}+\stackrel{}{k}^{}).$$
(45)
The Rayleigh-Jeans distribution corresponds to $`G(k)=T`$.

In the low-temperature phase, the long-wavelength part of the spectrum is suppressed by the Higgs effect,

$$G(k)=\frac{Tk^2}{k^2+m_\gamma ^2},$$
(46)
where $`m_\gamma `$ is the photon mass. These same considerations apply to superconductors, as well, and in that case $`m_\gamma `$ is the inverse penetration depth.

Following Ref. proc1 , one can calculate $`m_\gamma `$ as a function of temperature, and to leading order one finds,

$$m_\gamma ^2=e^2\varphi _0^2=\frac{e^2}{4}(T_c^2T^2).$$
(47)
Assuming again linear cooling, one finds

$$m_\gamma ^2=\frac{e^2T_c^2}{2}\frac{t}{\tau _Q}.$$
(48)
We will now assume that as the system is cooled through the transition, the scalar field remains close to equilibrium. This is not always true, and if it is not, then one would expect something like the Kibble mechanism. However, here we are interested in the opposite limit in which the gauge field falls out of equilibrium first. We also assume that the mean field result Eq. (48) is approximately valid.

According to Eq. (46), the amplitude $`G(k)`$ of a given mode decreases slowly at first, but when $`m_\gamma k`$, it starts to drop sharply toward zero. Using Eq. (48), we see that this happens at time

$$tt_k=\frac{2k^2}{e^2T_c^2}\tau _Q.$$
(49)
After this time, the equilibrium value of $`G(k)`$ decreases as

$$G_{\mathrm{eq}}(k)\frac{2k^2}{e^2T_c}\frac{\tau _Q}{t},$$
(50)
and we find

$$\frac{d\mathrm{ln}G_{\mathrm{eq}}(k)}{dt}=\frac{1}{t}.$$
(51)
This rate is faster at early times, and the fastest decay takes place at $`tt_k`$. Thus,

$$\frac{d\mathrm{ln}G_{\mathrm{eq}}(k)}{dt}|_{\mathrm{max}}\frac{1}{t_k}=\frac{e^2T_c^2}{2k^2\tau _Q}.$$
(52)
This is the rate at which $`G(k)`$ has to be able to decay for the mode to stay in equilibrium. We can see that long wavelengths (low $`k`$) have to decay faster.

On the other hand, one would expect on physical grounds that the dynamics of longer wavelengths is slower and that they would therefore actually decay shower than short wavelengths. It is therefore inevitable that some very long-wavelength modes are unable to reach the necessary rate in Eq. (52) and fall out of equilibrium instead.

We can be more specific if we assume that the dynamics is governed by the conductivity $`\sigma `$ and Ohm’s law is valid,

$$\stackrel{}{ȷ}=\sigma \stackrel{}{E}.$$
(53)
Maxwell’s equations then imply

$$\ddot{\stackrel{}{B}}=\stackrel{}{}\times \stackrel{}{}\times \stackrel{}{B}+\stackrel{}{}\times \stackrel{}{ȷ}=\stackrel{}{}\times \stackrel{}{}\times \stackrel{}{B}\sigma \dot{\stackrel{}{B}}.$$
(54)
In Fourier space (and ignoring vector indices), this becomes simply

$$\ddot{B}+\sigma \dot{B}+k^2B=0,$$
(55)
which is nothing but a damped harmonic oscillator. If $`k<\sigma /2`$, it has an exponential solution

$$B\mathrm{exp}\left(\frac{k^2}{\sigma }t\right).$$
(56)
For $`G(k)`$, this means that any deviation from the equilibrium value will decay no faster than at this rate,

$$\frac{d\mathrm{ln}G(k)}{dt}\frac{2k^2}{\sigma }.$$
(57)
Thus, we find that the critical value $`k_c`$, defined in such a way that modes with $`k<k_c`$ fall out of equilibrium, is given by

$$\frac{2k_c^2}{\sigma }=\frac{e^2T_c^2}{2k_c^2\tau _Q}k_c\left(\frac{e^2T_c^2\sigma }{4\tau _Q}\right)^{1/4}.$$
(58)
One can then roughly say that as $`m_\gamma \mathrm{}`$, modes with $`k>k_c`$ decay to zero, whereas those with $`k<k_c`$ retain their original amplitude. Thus, we would have the spectrum

$$G(k)=\{\begin{array}{cc}T,\hfill & k<k_c,\hfill \\ 0,\hfill & k>k_c.\hfill \end{array}$$
(59)
What this means is that the magnetic field has not disappeared completely. However, because the field cannot penetrate the broken, superconducting, phase, the magnetic field must be confined into flux tubes, i.e., vortices. Thus, we have seen that in the gauge theory, vortices are formed by the non-equilibrium dynamics of the gauge field. This phenomenon is knows as flux trapping IJMPA , because these vortices originate in the thermal fluctuations of the magnetic field, which do not have time to decay and are instead trapped in vortices.

In general, the winding number $`N_W`$ and the magnetic flux $`\mathrm{\Phi }=d^2\stackrel{}{S}\stackrel{}{B}`$ are related by the flux quantum $`\mathrm{\Phi }_0=2\pi /e`$, and for the case of a circle of radius $`R`$ in Eq. (39), we have

$$N_W(R)=\mathrm{\Phi }(R)/\mathrm{\Phi }_0e\mathrm{\Phi }(R).$$
(60)
On average, this is zero, but as in the global case, the variance is non-zero and can be calculated from $`G(k)`$,

$`N_W(R)^2`$ $``$ $`e^2\mathrm{\Phi }(R)^2`$ (61)
$`=`$ $`e^2{\displaystyle ^R}d^2xd^2yB_z(x)B_z(y)`$
$`=`$ $`e^2{\displaystyle ^R}d^2xd^2y{\displaystyle \frac{d^3k}{(2\pi )^3}e^{i\stackrel{}{k}(\stackrel{}{x}\stackrel{}{y})}G(k)}.`$
Using Eq. (59), this is

$$N_W(R)^2\{\begin{array}{cc}e^2T_cR,\hfill & \text{if }R>1/k_c,\hfill \\ e^2T_ck_c^3R^4,\hfill & \text{if }R<1/k_c.\hfill \end{array}$$
(62)
The scaling exponent $`\nu _D`$ is therefore

$$\nu _D=\{\begin{array}{cc}1,\hfill & \text{at short distances},\hfill \\ 1/4,\hfill & \text{at long distances}.\hfill \end{array}$$
(63)
Unfortunately, the long-distance value $`1/4`$ is equal to what the Kibble-Zurek mechanism gives, and cannot therefore be used to distinguish between the two theories. Because it depends on the dimensionality, the situation is different in two-dimensional models that have been studied numerically Hindmarsh:2000kd ; Stephens:2001fv . Note, however, that a two-dimensional film in three-dimensional space, the typical setup of actual experiments, behaves like a three-dimensional system in this sense Kibble:2003wt .

In three dimensions, the best way to distinguish between the theories is the short-distance exponent $`\nu =1`$. It means that the typical value of $`N_W(R)`$ grows as $`R^2`$, i.e., it is proportional to the area of the circle. This is only possible if all the vortices inside the circle have the same sign, or in other words, if there are clusters of equal-sign vortices (see Fig. 6). Eq. (62) tells us that the radius of these clusters is $`1/k_c`$, and the number of vortices in each is $`\sqrt{e^2T_c/k_c}`$. Consequently, the number density per unit cross sectional area is

$$n\frac{\sqrt{e^2T_c/k_c}}{1/k_c^2}\sqrt{e^2T_ck_c^3}\sqrt{e^2T}\left(\frac{e^2T_c^2\sigma }{\tau _Q}\right)^{3/8}.$$
(64)
Here we have only discussed Abelian gauge field theories, but similar arguments seem to apply to non-Abelian theories as well. In that case, one can also study the formation of magnetic monopoles Rajantie:2002dw , which may shed more light on the monopole problem in cosmology.

## XI Conclusions

We have seen that if a symmetry breaking phase transition takes place in a finite time, topological defects are formed. In a global theory, they are produced by the Kibble-Zurek mechanism, which leads to a characteristic negative correlation between vortices that can be used to identify the mechanism. In a gauge theory, defects are also produced by thermal fluctuations of the magnetic field that are too slow to decay. The clearest prediction of this flux trapping scenario is the formation of clusters of equal-sign vortices.

These theories are currently being tested in condensed matter experiments, in particular with superconductors. Although the calculations done in these paper refer to relativistic field theories, the qualitative conclusions would be the same in condensed matter systems.

###### Acknowledgements.

I would like to thank the organisers for a very successful workshop, and the ESF COSLAB programme and Churchill College for financial support.


text:
# 1 Introduction

## 1 Introduction

It was realized only recently that the nuclear modification factor defined as

$$R_{AA}^h(y,p_t;W)=\frac{\frac{d\sigma ^{AAh}}{dyd^2p_t}(y,p_t;W)}{N_{bin}\frac{d\sigma ^{pph}}{dyd^2p_t}(y,p_t;W)}$$
(1)
may be a very useful quantity which could shed more light on the dynamics of high-energy nuclear collisions . In the above definition $`y`$, $`p_t`$ and $`W`$ are hadron rapidity, transverse momentum, and nominal $`NN`$ center-of-mass energy, respectively. $`N_{bin}`$ is the number of binary nucleon-nucleon collisions usually calculated in the Glauber formalism. This brings in a small model dependence in the above definition. The index $`h`$ indicates the measured hadron.

Quite different values of nuclear modification factors have been observed at CERN SPS  and at the Brookhaven RHIC . It is speculated that the inequality

$$R_{AA}^h(RHIC)<R_{AA}^h(SPS)$$
(2)
is due to the creation at RHIC of a dense nuclear medium, a new state of matter, perhaps the long seeked quark-gluon plasma . Various nuclear effects like shadowing, parton multiple scattering or final state quenching  modify the dependence of $`R_{AA}`$ on the transverse momentum and collision energy in quite different ways. In Ref. the data from CERN SPS were interpreted within the picture of Glauber multiple scattering of initial partons (see ). It is not easy to make a priori predictions based on this effect as one has to model parton-nucleon interaction. The proposed schematic models are flexible enough so that a rough description of the data is possible.

In a partonic approach improved for including internal transverse momenta , the average transverse momentum of partons in nuclear collisions is increased compared to the one in the elementary collisions

$$p_t^2_{pA}=p_t^2_{pp}+\mathrm{\Delta }p_t^2(A,b,\mathrm{}).$$
(3)
The correction term used in the phenomenological analyses is typically a fraction of GeV, i.e. of the order of the nuclear Fermi momentum.

At intermediate energies the inclusion of momentum distributions in nuclei quantifying the nucleon motion is crucial for the production of kaons , eta mesons  and antiprotons  in p+A collisions. It is well known that at subthreshold and near-to-threshold collisions the Fermi motion leads to a huge enhancement of hadron production cross sections based on the first collision model. At the lower energies the inclusion of the off-shell effects incorporated through spectral functions (see for example ) can be numerically important.

At large energies one is tempted to completely ignore the effect of the Fermi motion of nucleons in the colliding nuclei. According to our knowledge, however, no quantitative estimates have been presented so far in the literature.

In the present Letter we wish to show that the Fermi motion effect on $`R_{AA}(p_t)`$ is not negligible even at high energies (W $`>`$ 15 GeV) and must be included into consideration, e.g. in the context of the discussion on the creation of quark-gluon plasma.

## 2 Folding model

In order to demonstrate our main point, in the following we shall neglect all nuclear effects except the nucleon Fermi motion. We assume that the collision of nuclei $`A_1(Z_1,N_1)`$ and $`A_2(Z_2,N_2)`$ can be considered as the collision of two ensambles of nucleons. In general, the production cross section of a hadron $`h`$ depends on the type of the collision ($`\alpha ,\beta `$ = pp, pn, np or nn). We assume that at these high energies all off-shell effects can be neglected. In this somewhat simplified picture the cross section per elementary collision 111This corresponds to $`\frac{\frac{d\sigma ^{AAh}}{dyd^2P_t}}{N_{bin}}`$ in the definition of the nuclear modification factor (1). is

$$\begin{array}{cc}& \frac{d\sigma ^{AAh}}{dyd^2p_t}(y,p_t;W)d^3p_1d^3p_2\underset{\alpha ,\beta =p,n}{}\frac{N_{\alpha ,\beta }}{N_{tot}}f_{\alpha /1}(p_1)f_{\beta /2}(p_2)\hfill \\ & 𝑑y^{}d^2p_t^{}\delta (\stackrel{~}{y}y)\delta ^2(\stackrel{}{\stackrel{~}{p}}_t\stackrel{}{p}_t)\frac{d\sigma ^{\alpha \beta h}}{dy^{}d^2p_t^{}}(y^{},p_t^{};W^{}).\hfill \end{array}$$
(4)
In the above equation $`y^{}`$, $`p_t^{}`$, $`W^{}`$ denote rapidity, transverse momentum and energy in the nucleon-nucleon center-of-mass system, respectively. The quantities $`\stackrel{~}{y}`$ and $`\stackrel{}{\stackrel{~}{p}_t}`$ (overall CM system) are obtained from $`\stackrel{}{p}_1`$, $`\stackrel{}{p}_2`$, $`y^{}`$, $`\stackrel{}{p}_t^{}`$, etc., via appropriate Lorentz transformations and rotations of the frame of reference. In Eq.(4) $`p_1`$ is the nucleon momentum in the rest frame of the nucleus $`A_1`$, while $`p_2`$ is the nucleon momentum in the rest frame of the nucleus $`A_2`$. Eq.(4) is valid provided the standard normalization of momentum distributions is used

$$f_{\alpha /i}(p_i)d^3p_i=1.$$
(5)
In the master formula (4), $`N_{\alpha ,\beta }`$ is the effective number of the collisions between nucleons $`\alpha `$ and $`\beta `$. Obviously $`N_{tot}=_{\alpha ,\beta }N_{\alpha ,\beta }`$. In the most naive approach $`N_{\alpha ,\beta }/N_{tot}=(N_{\alpha /1}N_{\beta /2})/(A_1A_2)`$.

In our intentionally simplified approach we consider the nucleus-nucleus collision as a collision of two bunches of nucleons. For simplicity we make no distinction between proton-proton, proton-neutron, neutron-proton and neutron-neutron collisions 222This means that our present approach is valid only for the cases when: $`\sigma _{pph}=\sigma _{pnh}=\sigma _{nph}=\sigma _{nnh}`$. This is roughly true for example for h= $`\pi ^0`$, $`\pi ^++\pi ^{}`$, $`\eta `$.. In addition, we shall make simplifications of kinematics relevant at high energies, when collective velocities of both nuclei are much larger than the Fermi velocity, i.e. $`v_{A_1/CM},v_{A_2/CM}v_F`$. In this case $`\stackrel{}{v}_{1/NN},\stackrel{}{v}_{2/NN}\widehat{z}\stackrel{}{v}_{A_1/CM},\stackrel{}{v}_{A_2/CM}`$ and $`y^{}y`$. Then, the invariant distribution (per elementary collision) of hadron $`h`$ in nuclear collisions is given as

$$\frac{d\sigma ^{AAh}}{dyd^2p_t}(y,p_t;W)d^3p_1d^3p_2f_1(p_1)f_2(p_2)\frac{d\sigma ^{NNh}}{dy^{}d^2p_t^{}}(y^{},p_t^{};W^{}).$$
(6)
In the formula above, $`W^{}`$ denotes the energy available in the elementary nucleon-nucleon collision. In order to calculate $`W^{}`$ one needs to make the appropriate Lorentz transformations from four momenta in the rest frame of nuclei $`(E_{i/A_i},\stackrel{}{p}_{i/A_i})`$ to four-momenta in the overall center of mass of the colliding nuclei $`(E_{i/CM},\stackrel{}{p}_{i/CM})`$. Then the nucleon-nucleon c.m. energy is calculated as

$$W^{}=\sqrt{\left(E_{1/CM}+E_{2/CM}\right)^2\left(\stackrel{}{p}_{1/CM}+\stackrel{}{p}_{2/CM}\right)^2}.$$
(7)
In the equation above the on-shell relation $`E_{i/CM}=\sqrt{p_{i/CM}^2+m_N^2}`$ ($`i`$=1,2) is used. Due to plateau at midrapidities, at $`y`$ 0 it seems reasonable to make an extra approximation $`y^{}=y`$.

Although the master formula (6) is very simple, its practical application requires knowledge of nucleon momentum distributions ($`f_1`$ and $`f_2`$) in nuclei as well as the rapidity, transverse momentum and energy dependence of the elementary production mechanism. It is needless to say that realistic calculation of both nucleon momentum distribution in nuclei and differential cross section for elementary collisions is rather difficult. While the first requires nuclear many-body approach to the nucleus structure, the latter involves nonperturbative QCD effects which are not fully understood. The exact nuclear calculation based on modern nucleon-nucleon interactions are available only for the deuteron, $`{}_{}{}^{3}H`$, $`{}_{}{}^{3}He`$, and quite recently – the four-nucleon system $`{}_{}{}^{4}He`$ . Similarly, our understanding of the soft sector (small transverse momenta) of particle production within QCD is still not complete .

## 3 Estimate of the effect

Before we go to the presentation of predictions based on Eq.(6), we shall shortly describe the ingredients which enter the formula. In the Fermi gas model

$$f(p)=C_{FG}\theta (p_Fp),$$
(8)
where $`p_F`$ is Fermi momentum and $`\theta `$ is the standard step function. The normalization constant $`C_{FG}`$ can be obtained from the normalization condition (5). In numerical calculations we take $`p_F`$ = 1.37 fm-1. In the present paper we shall also consider a simple Gaussian parametrization of the momentum distribution of nucleons in $`{}_{}{}^{12}C`$ nucleus from Ref.. While the first seems more adequate for heavy nuclei, the latter may be representative for light nuclei.

The elementary cross sections for the production of $`\pi ^\pm ,\pi ^0,K^\pm `$ in the proton-proton or proton-antiproton collisions are know experimentally in the broad range of energies. As mentioned above, in the present analysis we shall limit ourselves to the production of $`\pi ^0`$ only. For our purpose a simple parametrization of the available world data would be useful. Such a parametrization exists, however, only for charged particles in the limited range of energies . In the parton picture, assuming isospin symmetry, the invariant cross section for the $`pp\pi ^0X`$ reaction can be expressed as

$$\frac{d\sigma ^{pp\pi ^0}}{dyd^2p_t}(y,p_t;W)=\frac{1}{2}\left(\frac{d\sigma ^{pp\pi ^+}}{dyd^2p_t}(y,p_t;W)+\frac{d\sigma ^{pp\pi ^{}}}{dyd^2p_t}(y,p_t;W)\right).$$
(9)
The effects which violate the above relation are expected to be sufficiently small. In the present calculation we use a parametrization of the ISR data from Ref..

In general, there are three effects contained in the original convolution formula (4) which may result in Fermi motion modifying of $`R_{AA}^h`$. The first effect is connected to the fact that $`p_t^{}p_t`$. The second effect is related to the fact that $`y^{}y`$. Because the general Lorentz transformation from the NN system to the overall nucleus-nucleus CM system includes also rotations of the frame axes, the formulae are rather cumbersome and not easy to reverse. Therefore these two effects cannot be included exactly in a simple way. The Monte Carlo method seems more suitable for this purpose. A Monte Carlo calculation, to be discussed in details elsewhere, shows that while the first effect is rather small, the second effect is almost negligible.

The main (third) effect is connected to the fact that the true c.m. energy in the elementary nucleon-nucleon collision is not equal to the nominal nucleon-nucleon energy $`W^{}W`$, even when averaged over Fermi distributions. As an example, in Fig.1 we show two distributions of the true nucleon-nucleon center of mass energy $`f(W^{})`$ 333$`f(W^{})𝑑W^{}`$ = 1. for two nominal energies $`W`$ = 20 GeV (left panel) and $`W`$ = 50 GeV (right panel). We observe a broadening of $`W^{}`$ distribution with raising nominal energy, although the relative dispersion $`\mathrm{\Delta }W^{}/W`$ stays roughly constant. The invariant elementary cross sections for particle production rise across the region over which $`f(W^{})`$ is spread (see Fig.2) leading to the value of $`W^{}_\sigma `$ 444The index $`\sigma `$ indicates averaging of $`W^{}`$ with the energy-dependent $`NNh`$ cross section. larger than $`W`$. This effect is stronger for larger $`p_t^{}`$ and smaller energies $`W^{}`$.

In Fig.3 we show $`R_{AA}^{\pi ^0}`$ as a function of the $`\pi ^0`$ transverse momentum for the $`\pi ^0`$ rapidity y=0 and two different Fermi momentum distributions:
(a) Fermi gas model distribution given by Eq.(8) (dashed) and
(b) the Gaussian parametrization from Ref. (solid).
The experimental data points were obtained by dividing the nuclear data from Ref. by the cross section for the elementary collision obtained with the help of Eq.(9) and the parametrization from .

The Fermi motion effect calculated from our master formula (6) is strongly energy dependent, which is demonstrated in Fig.4. Generally, the larger energy the smaller the effect. Surprisingly, the effect of the Fermi motion stays present even up to W = 50 GeV.

In the present analysis we have estimated the effect of Fermi motion assuming that nucleons can be treated as free on-mass-shell particles. Under this assumption the light-cone fraction carried by a nucleon, $`\xi `$, defined so that 0$`<\xi <`$A, is $`\xi =(\sqrt{m_N^2+p^2}+p_3)/(m_A/A)`$. As a result the average $`<\xi >`$ is larger than 1 and the momentum sum rule is violated - nucleons carry momenta larger than the total nucleus momentum. Effectively the energy in NN collisions is overestimated. This is connected with the effect of nucleon binding in nuclei.

The problem can be even better seen in the nucleus rest frame, where the total nucleon-gas energy is larger than the nucleus mass

$$A\sqrt{p^2+m_N^2}f(p)d^3p>M_A,$$
(10)
i.e. the gas of nucleons has too large energy. This effect can be cured introducing effective nucleon mass $`m_N^{eff}`$:

$$A\sqrt{p^2+(m_N^{eff})^2}f(p)d^3p=M_A.$$
(11)
This corresponds to the reduction of the nucleon mass by 10-20%.

Can the binding effect change our conclusions? The binding effect in nuclei is in general a complicated multiparticle problem and was solved in the exact way essentially only for three- and four-nucleon systems. The binding effect in heavy nuclei can be approximately included by introducing a notion of nucleon effective mass in the relation $`E_{N/A}=\sqrt{p_{N/A}^2+(m_N^{eff})^2}`$. Since the binding effect is only a correction to our main effect of Fermi motion this approach seems justified. For illustration in Fig.5 we show the result for nominal energy W = 20 GeV with nucleon mass modified by $`\pm `$ 150 MeV. This figure shows that the binding effect is not only qualitatively but also quantitatively a correction to our effect. This correction does not change the main trend of $`R_{AA}(p_t)`$. A more refined calculation is not our intention and clearly goes beyond the scope of the present paper.

How does the effect of Fermi motion depend on the type of particle measured? In Fig.6 we present $`R_{AA}`$ for $`\pi ^++\pi ^{}`$ (solid), $`K^++K^{}`$ (dashed) and $`p+\overline{p}`$ (dotted) at the center-of-mass energy $`W`$ = 20 GeV and $`y`$ = 0. In this calculation Fermi gas model was used. The results for pions and kaons are similar, while the result for protons+antiprotons are somewhat different. The latter does not show the enhancement at large transverse momenta observed for pions and kaons. This is due to the parametrization for protons and antiprotons  which is practically energy independent. Perhaps the difference for protons+antiprotons and pions/kaons can be used to identify experimentally the Fermi motion effect.

In the present letter we have concentrated only on the effect of Fermi motion on nuclear modification factor. The Fermi motion in colliding nuclei is probably important in many other fixed-target (low energy) processes. For example global partonic QCD analyses use Drell-Yan data based on rather low-energy proton-nucleus collisions than elementary ones. The standard QCD approach is known to have problems in describing the transverse momentum distribution of the lepton pair. In such analyses the Fermi motion of nucleons in nuclei is completely ignored. We expect that the effect of the Fermi motion should be dependent on the invariant mass of the two leptons. In general, the larger invariant mass, the larger Fermi motion effect may be expected. From the experience of the present paper one may expect large effects close to kinematical limits. Therefore even for the production of heavy gauge bosons (W,Z) in nucleus-nucleus collisions at RHIC one may expect sizeable effect due to the Fermi motion. These effects were never discussed in the literature and await a quantitative analysis.

Finally we wish to mention another interesting, not fully explained, effect related to the CERES collaboration measurement of the elliptic flow parameter $`v_2`$ . The existing hydrodynamical or thermodynamical models have serious problems in explaining a saturation of $`v_2`$ observed at larger transverse momenta (see e.g.). It is quite probably that the Fermi motion effect may play important role in the full understanding of this effect. However, this effect is very subtle and requires the inclusion of many effects in addition to the nucleon Fermi motion. This clearly goes beyond the scope of the present, intentionally simplified, analysis.

## 4 Conclusions

For the first time in the literature, in the present Letter we have analysed the role of the Fermi motion in nuclei in modyfing the so-called nuclear modification factor $`R_{AA}`$. A large effect is found even at SPS energies. The details depend on the model of the Fermi motion as well as on the cross sections for elementary collisions. The latter are very difficult to calculate reliably, especially at low transverse momenta. Therefore we have made use of a parametrization of existing experimental data.

At SPS energies the Fermi motion causes a sizeable enhancement at $`p_t>`$ 4 GeV. At RHIC energies the enhancement almost dissapears.

This quantitative difference of the role of the Fermi motion at low and high energies is very important in the context of a quantitative understanding of jet quenching in dense nuclear media like quark-gluon plasma. At present the latter cannot be calculated from first principles and the model parameters are extracted by comparison to the experimental data.

Even when alone, the Fermi motion leads to the Cronin effect. Different mechanisms like parton multiple scattering , initial transverse momentum broadening in nuclei  or perturbative saturation of gluon distributions  can contribute to the Cronin effect too. Our analysis shows that the Fermi motion effect is very important and must be included into the analysis before reliable conclusions about other nuclear effects may be drawn from such phenomenological studies. We find $`R_{AA}^{\pi ^0}(RHIC)<R_{AA}^{\pi ^0}(SPS)`$ which makes the interpretation of the experimental data even more complicated. For example at the SPS energy W=17.4 GeV the inclusion of the Fermi motion would improve the agreement of the calculation in Ref. with the WA98 collaboration experimental data (see their Fig.3). However, the inclusion of all effects like shadowing, quenching and Fermi motion in one consistent framework seems rather difficult at present.

The Fermi motion effect seems to be one of a few different competing mechanisms modyfing $`R_{AA}`$. Can one isolate experimentally the Fermi motion effect better? We think that the role of Fermi motion can be pinned down better in the collisions of light nuclei, where other effects are suppressed and the nucleon momentum distribution can be calculated from realistic models of nucleon-nucleon interactions. The examples are $`dd`$ or $`\alpha \alpha `$ scattering measured already in the past (see for example ). Better precision data at SPS energies would help to identify the effect of the Fermi motion via a comparison with realistic theoretical calculations . Also very precise data on elementary collisions would be helpful when calculating the modification of $`R_{AA}`$ precisely.

Acknowledgments We are indebted to Klaus Reygers for providing us with the WA98 collaboration experimental data for $`\pi ^0`$ production in $`{}_{}{}^{208}Pb+^{208}Pb`$ collisions and Zbigniew Rudy for providing us with a parametrization for the nucleon momentum distribution in $`{}_{}{}^{12}C`$. One of us (A.S.) is indebted to Piotr Bożek and Wojciech Florkowski for an interesting discussion and Piotr Żenczykowski for his careful reading the manuscript. We are indebted to Kolya Nikolaev and Andrzej Rybicki for a discussion which helped in finding a mistake in our code and therefore helped in correcting a first version of our letter.