text:
# Manipulating the Quantum State of an Electrical Circuit

(Dec. 21, 2001)

## Abstract

We have designed and operated a superconducting tunnel junction circuit that behaves as a two-level atom: the “quantronium”. An arbitrary evolution of its quantum state can be programmed with a series of microwave pulses, and a projective measurement of the state can be performed by a pulsed readout sub-circuit. The measured quality factor of quantum coherence $`Q_\phi 25000`$ is sufficiently high that a solid-state quantum processor based on this type of circuit can be envisioned.

Can we build machines that actively exploit the fundamental properties of quantum mechanics, such as the superposition principle or the existence of entangled states? Applications such as the transistor or the laser, often quoted as developments based on quantum mechanics, do not actually answer this question. Quantum mechanics enters into these devices only at the level of material properties but their state variables such as voltages and currents remain classical. Proposals for true quantum machines emerged in the last decades of the 20th century and are now being actively explored: quantum computers QC , quantum cryptography communication systems QS2  and detectors operating below the standard quantum limit QA3 . The major difficulty facing the engineer of a quantum machine is decoherence decoh2 . If a degree of freedom needs to be manipulated externally, as in the writing of information, its quantum coherence usually becomes very fragile. Although schemes that actively fight decoherence have recently been proposed QEC ; QEC2 , they need very coherent quantum systems to start with. The quality of coherence for a two-level system can be quantitatively described by the quality factor of quantum coherence $`Q_\phi =`$ $`\pi \nu _{01}T_\phi `$ where $`\nu _{01}`$ is its transition frequency and $`T_\phi `$ is the coherence time of a superposition of the states. It is generally accepted that for active decoherence compensation mechanisms, $`Q_\phi `$’s larger than $`10^4\nu _{01}t_{\mathrm{op}}`$ are necessary, $`t_{\mathrm{op}}`$ being the duration of an elementary operation minQ .

Among all the practical realizations of quantum machines, those involving integrated electrical circuits are particularly attractive. However, unlike the electric dipoles of isolated atoms or ions, the state variables of a circuit like voltages and currents usually undergo rapid quantum decoherence because they are strongly coupled to an environment with a large number of uncontrolled degrees of freedom Schoen . Nevertheless, superconducting tunnel junction circuits MQT2 ; Nakamura ; Mooij ; Lukens ; Han  have displayed $`Q_\phi `$’s up to several hundred Martinis  and temporal coherent evolution of the quantum state has been observed on the nanosecond time scale Nakamura ; Nakamura3  in the case of the single Cooper pair box quantro . We report here a new circuit built around the Cooper pair box with $`Q_\phi `$ in excess of $`10^4`$, whose main feature is the separation of the write and readout ports design ; Zorin . This circuit, which behaves as a tunable artificial atom, has been nicknamed a “quantronium”. The basic Cooper pair box consists of a low capacitance superconducting electrode, the “island”, connected to a superconducting reservoir by a Josephson tunnel junction with capacitance $`C_j`$ and Josephson energy $`E_J`$. The junction is biased by a voltage source $`U`$ in series with a gate capacitance $`C_g`$. In addition to $`E_J`$ the box has a second energy scale, the Cooper pair Coulomb energy $`E_{CP}=(2e)^2/2\left(C_g+C_j\right)`$. When the temperature $`T`$ and the superconducting gap $`\mathrm{\Delta }`$ satisfy $`k_BT\mathrm{\Delta }/\mathrm{ln}𝒩`$ and $`E_{CP}\mathrm{\Delta }`$, where $`𝒩`$ is the total number of paired electrons in the island, the number of excess electrons is even Tuominen ; parity . The Hamiltonian of the box is then

$$\widehat{H}=E_{CP}\left(\widehat{N}N_g\right)^2E_J\mathrm{cos}\widehat{\theta }\text{,}$$
(1)
where $`N_g=C_gU/2e`$ is the dimensionless gate charge and $`\widehat{\theta }`$ the phase of the superconducting order parameter in the island, conjugate to the number $`\widehat{N}`$ of excess Cooper pairs in it quantro .

In our experiment, $`E_JE_{CP}`$ and neither $`\widehat{N}`$ nor $`\widehat{\theta }`$ is a good quantum number. The box thus has discrete quantum states that are quantum superpositions of several charge states with different $`N`$. Because the system is sufficiently non-harmonic, the ground $`|0`$ and first excited $`|1`$ energy eigenstates form a two-level system. This system corresponds to an effective spin one-half $`\stackrel{}{s}`$ whose Zeeman energy $`h\nu _{01}`$ goes to a minimal value close to$`E_J`$ when $`N_g=1/2`$. At this particular bias point both states $`|0`$($`s_z=+1/2`$) and $`|1`$ ($`s_z=1/2`$) have the same average charge $`\widehat{N}=1/2,`$ and consequently the system is immune to first order fluctuations of the gate charge. Manipulation of the quantum state is performed by applying microwave pulses $`u(t)`$ with frequency $`\nu \nu _{01}`$ to the gate, and any superposition $`|\mathrm{\Psi }=\alpha |0+\beta |1`$ can be prepared.

A novel type of readout has been implemented in this work. The single junction of the basic Cooper pair box has been split into two nominally identical junctions in order to form a superconducting loop (Fig. 1). The Josephson energy $`E_J`$ in Eq. 1 becomes $`E_J\mathrm{cos}(\widehat{\delta }/2)`$ Likharev , where $`\widehat{\delta }`$ is an additional degree of freedom, the superconducting phase difference across the series combination of the two junctions Averin . The two states are discriminated not through the charge $`\widehat{N}`$ on the island Nakamura ; Delsing , but through the supercurrent in the loop $`\widehat{I}=(2e/\mathrm{})\widehat{H}/\widehat{\delta }`$. This is achieved by entangling $`\stackrel{}{s}`$ with the phase $`\widehat{\gamma }`$ of a large Josephson junction with Josephson energy $`E_{J0}20E_J`$, inserted in the loop design ; Buisson . The phases are related by $`\widehat{\delta }=\widehat{\gamma }+`$ $`\varphi `$, where $`\varphi =2e\mathrm{\Phi }/\mathrm{}`$, $`\mathrm{\Phi }`$ being the external flux imposed through the loop. The junction is shunted by a capacitor $`C`$ to reduce phase fluctuations. A trapezoidal readout pulse $`I_b(t)`$ with a peak value slightly below the critical current $`I_0=2eE_{J0}/\mathrm{}`$ is applied to the parallel combination of the large junction and the small junctions (Fig. 1C). When starting from $`\widehat{\delta }0`$, the phases $`\widehat{\gamma }`$ and $`\widehat{\delta }`$ grow during the current pulse, and consequently a $`\stackrel{}{s}`$-dependent supercurrent develops in the loop. This current adds to the bias-current in the large junction, and by precisely adjusting the amplitude and duration of the $`I_b(t)`$ pulse, the large junction switches during the pulse to a finitevoltage state with a large probability $`p_1`$ for state $`|1`$ and with a small probability $`p_0`$ for state $`|0`$design . This readout scheme is similar to the spin readout of Ag atoms in a Stern and Gerlach apparatus, in which the spin is entangled with the atom position. For the parameters of the experiment, the efficiency of this projective measurement should be $`\eta =p_1p_0=0.95`$ for optimum readout conditions. The readout is also designed so as to minimize the $`|1|0`$ relaxation rate using a Wheatstone-bridge-like symmetry. Large ratios $`E_{J0}/E_J`$ and $`C/C_j`$ provide further protection from the environment. Just as the system is immune to charge noise at $`N_g=1/2`$, it is immune to flux and bias current noise at $`\varphi =0`$ and $`I_b=0,`$ where $`\widehat{I}=0`$. The preparation of the quantum state and its manipulation are therefore performed at this optimal working point.

A quantronium sample is shown in Fig. 1B. It was fabricated with standard e-beam lithography and aluminum evaporation. The sample was cooled down to 15 mK in a dilution refrigerator. The switching of the large junction plasma  to the finite voltage state is detected by measuring the voltage across it with a room temperature preamplifier followed by a discriminator. By repeating the experiment, the switching probability, and hence the occupation probabilities of the $`|0`$ and $`|1`$ states, can be determined.

The readout part of the circuit was tested by measuring the switching probability $`p`$ at thermal equilibrium as a function of the pulse height $`I_p`$, for a readout pulse duration of $`\tau _r=100\mathrm{ns}`$. The discrimination between the estimated currents for the $`|0`$ and $`|1`$ states was found to have an efficiency of $`\eta =0.6`$, which is lower than the expected $`\eta =0.95`$. Measurements of the switching probability as a function of temperature and repetition rate indicate that the discrepancy between the theoretical and experimental readout efficiency could be due to an incomplete thermalization of our last filtering stage in the bias current line.

Spectroscopic measurements of $`\nu _{01}`$were performed by applying to the gate a weak continuous microwave irradiation suppressed just before the readout pulse. The variations of the switching probability as a function of the irradiation frequency display a resonance whose center frequency evolves with dc gate voltage and flux as the Hamiltonian predicts, reaching $`\nu _{01}16.5\mathrm{GHz}`$ at the optimal working point (Fig. 2). The small discrepancy between theoretical and experimental value of the transition frequency at non-zero magnetic flux is attributed to flux penetration in the small junctions not taken into account in the model. These spectroscopic data have been used to precisely determine the relevant circuit parameters, $`E_J=0.865k_B\mathrm{K}`$ and $`E_J/E_{CP}=1.27`$. At the optimal working point, the linewidth was found to be minimal with a $`0.8`$MHz full width at half-maximum (FWHM). When varying the delay between the end of the irradiation and the readout pulse, the resonance peak height decays with a time constant $`T_1=1.8\mu \mathrm{s}`$. Supposing that the energy relaxation of the system is only due to the bias circuitry, a calculation similar to that in relax  predicts that $`T_110\mu \mathrm{s}`$ for a crude discrete element model. This result shows that no detrimental sources of dissipation have been seriously overlooked in our circuit design.

Controlled rotations of $`\stackrel{}{s}`$ around an axis $`x`$ perpendicular to the quantization axis $`z`$ have been performed. Before readout, a single pulse at the transition frequency with variable amplitude $`U_{\mu \mathrm{w}}`$ and duration $`\tau `$ was applied. The resulting change in switching probability is an oscillatory function of the product $`U_{\mu \mathrm{w}}\tau `$ (Fig. 3A), which is in agreement with the theory of Rabi oscillations Rabi , proving that the resonance indeed arises from a two-level system. The proportionality ratio between the Rabi period and $`U_{\mu \mathrm{w}}\tau `$ was used to calibrate microwave pulses for the application of controlled rotations of $`\stackrel{}{s}`$.

Rabi oscillations correspond to a driven coherent evolution but do not give direct access to the intrinsic coherence time $`T_\phi `$ during a free evolution of $`\stackrel{}{s}`$ . This $`T_\phi `$ was obtained by performing a Ramsey fringe experiment Ramsey  on which atomic clocks are based. One applies to the gate two phase coherent microwave pulses each corresponding to a $`\pi /2`$ rotation around $`x`$ noteangle  and separated by a delay $`\mathrm{\Delta }t`$ during which the spin precesses freely around $`z`$. For a given detuning of the microwave frequency, the observed decaying oscillations of the switching probability as a function of $`\mathrm{\Delta }t`$ (Fig. 3B) correspond to the “beating” of the spin precession with the external microwave field noteRamsey . The oscillation period agrees exactly with the inverse of the detuning, allowing a measurement of the transition frequency with a relative accuracy of $`6\times 10^6`$. The envelope of the oscillations yields the decoherence time $`T_\phi 0.50\mu \mathrm{s}`$. Given the transition period $`1/\nu _{01}60\mathrm{ps}`$, this means that $`\stackrel{}{s}`$ can perform on average 8000 coherent free precession turns.

In all the time domain experiments on the quantronium, the oscillation period of the switching probability agrees closely with theory, which proves controlled manipulation of $`\stackrel{}{s}`$. However, the amplitude of the oscillations is smaller than expected by a factor of 3 to 4. This loss of contrast is likely to be due to a relaxation of the level population during the measurement itself.

In order to understand what limits the coherence time of the circuit, measurements of the linewidth $`\mathrm{\Delta }\nu _{01}`$ of the resonant peak as a function of $`U`$ and $`\mathrm{\Phi }`$ have been performed. The linewidth increases linearly when departing from the optimal point $`(N_g=1/2,\varphi =0,I_b=0)`$. This dependence is well accounted for by charge and phase noises with root mean square deviations $`\mathrm{\Delta }N_g=0.004`$ and $`\mathrm{\Delta }(\delta /2\pi )=0.002`$ during the time needed to record the resonance. The residual linewidth at the optimal working point is well-explained by the second order contribution of these noises. The amplitude of the charge noise is in agreement with measurements of $`1/f`$ charge noise PTB , and its effect could be minimized by increasing the $`E_J/E_{CP}`$ ratio. The amplitude of the flux noise is unusually large Wellstood  and should be significantly reduced by improved magnetic shielding. An improvement of $`Q_\phi `$ by an order of magnitude thus seems possible. Experiments on quantum gates based on the controlled entanglement of several capacitively coupledquantronium circuits could be already performed with the level of quantum coherence achieved in the present experiment.

###### Acknowledgements.

The indispensable technical work of P. Orfila is gratefully acknowledged. This work has greatly benefited from direct inputs from J. M. Martinis and Y. Nakamura. The authors acknowledge discussions with P. Delsing, G. Falci, D. Haviland, H. Mooij, R. Schoelkopf, G. Schön and G. Wendin. Partly supported by the European Union through contract IST-10673 SQUBIT and the Conseil Général de l’Essonne through the EQUM project.


format: text

text:
# Exotic fluids and crystals of soft polymeric colloids

## 1 Introduction

One major advantage of soft matter systems in comparison to their atomic counterparts is the ability one has to engineer the constituent particles at the molecular level. In this way, an enormous variety in architectures can be achieved, leading to a corresponding richness in the structural and phase behaviour of such systems. Polymers play a prominent example within this class of materials. They come in a variety of forms, such as, e.g., as linear chains, branched, star-shaped, dendritic, copolymer, as well as functions, e.g., steric stabilisers, additives, depletants etc. Additional flexibility arises from the possibility of influencing the structural and phase behaviour of polymer solutions by changing the solvent quality.

From the theoretical point of view, the task of bridging the gap between the microscopic and macroscopic length scales of polymeric systems is a formidable one. Indeed, the constituent macromolecules may contain thousands or even millions of atoms, interconnected to one another in complicated ways. The program of starting with the individual interactions between monomers and proceeding to the calculation of the free energy of the system rapidly becomes intractable. However, considerable progress can be made if one invokes a “coarse-graining” procedure, a point of view that has been proved very fruitful in many areas of research in condensed matter physics. Here we follow a two-step procedure. First, instead of attempting to carry out the statistical trace over all the individual monomers in one step, a certain generalised coordinate of the macromolecule as an effective point particle is invoked. Examples include the centre of mass of the polymer or some suitably selected monomer, as we explain below. All monomers belonging to the macromolecules are then traced out for a given, fixed configuration of the effective coordinates, which defines an effective interaction between these coordinates . Once this is achieved, the second step consists of viewing the the macroscopic system as collection of point particles interacting by means of the effective interaction. Now all known tools from the theory of atomic and molecular fluids can be employed to derive structural and thermodynamic quantities for the system under consideration.

In the recent years, the program sketched above has been carried out with success for various polymeric systems. It has been found that the effective interactions obtained belong to a new class of ultrasoft potentials which have very unusual properties when compared with the hard-sphere (HS) system, the prototype of atomic liquids. In this work, we first present a concise review of these novel properties for the fluid phases and then some new results regarding the rich variety of crystal structures that can be stabilised by ultrasoft potentials. The rest of the paper is organised as follows. In section 2 we give the general definition of the effective interaction and present specific examples for a number of systems that have been worked out recently. In section 3 we show that the systems described by ultrasoft interactions are well described by a simple mean-field picture in the fluid phase for a wide range of thermodynamic conditions. To further delineate these systems we define two categories of mean-field fluids, the strong and the weak ones. The former are characterised by a direct correlation function $`c(r)`$ that satisfies to excellent accuracy the condition $`c(r)=\beta v(r)`$ in broad regions of the thermodynamics space, where $`v(r)`$ is the interparticle pair potential. The latter only satisfy an approximate mean-field picture for their thermodynamic properties, and not for their structure. In section 4 we demonstrate that strong mean-field fluids can be further divided into two categories: those displaying reentrant melting and those displaying a cascade of clustering transitions, the criterion being set by the positivity of the Fourier transform of the effective interaction. In section 5 we discuss the richness of the phase diagrams of weak mean-field fluids, taking the case of star polymers as a concrete example. Reentrant melting as well as a wealth of structural phase transitions and exotic crystal phases are all shown to be stemming from the ultrasoftness of the effective interactions. Finally, in section 6 we summarise and conclude.

## 2 Effective interactions between polymeric macromolecules

The effective interaction between flexible, fluctuating aggregates can be given a precise statistical-mechanical definition. Let us consider a solution containing $`M`$ polymeric macromolecules, each one of them composed of $`N`$ monomers. The total number of particles in the system is $`=M\times N`$. One starts from the full Hamiltonian $``$ of the problem, assumed to be known. Then, out of the $``$ particles in the problem (in our case all monomers), one selects the $`M`$ ones that are to be considered as “effective particles” and holds them fixed in some prescribed configuration $`\{𝐑_1,𝐑_2,\mathrm{},𝐑_M\}`$, where $`𝐑_i`$ is the position of the $`i`$-th effective particle. Afterwards, the $`M`$ remaining particles are canonically traced out and the result of this integration is a constrained partition function $`Q(𝐑_1,𝐑_2,\mathrm{},𝐑_M)`$. The effective Hamiltonian $`_{\mathrm{eff}}`$ is defined as:

$$\mathrm{exp}\left(\beta _{\mathrm{eff}}\right)=Q(𝐑_1,𝐑_2,\mathrm{},𝐑_M),$$
(1)
where $`\beta =(k_BT)^1`$, with the absolute temperature $`T`$ and Boltzmann’s constant $`k_B`$. It can be shown  that such an effective Hamiltonian has two useful properties: it preserves the overall thermodynamics of the system and it guarantees that the correlation functions of any order between any of the $`M`$ remaining particles remain invariant, regardless of whether the expectation values are calculated with the original Hamiltonian $``$ or with the effective one $`_{\mathrm{eff}}`$. Though the procedure of tracing out the $`M`$ degrees of freedom necessarily generates interactions of all orders between the $`M`$ particles  in many cases it is sufficient to truncate those at the pair level, introducing thereby the pair-potential approximation. The great advantage of employing this point of view is that, in comparison with the original problem, the numbers of particles, $``$, has been reduced by a factor $`N`$. In addition, whereas in the original problem the pair interactions between the monomers are quite complicated, due to the need of taking into account the connectivity and architecture of the molecule, in the effective description the pair potential is spherically symmetric, depending only on the magnitude of the separation vector between the two effective coordinates. A new picture of the polymers emerges thereby in which the latter can be seen as ultrasoft colloids having a dimension of the order of their radius of gyration $`R_g`$. This sets at the same time the characteristic length scale of the effective interaction between them. When the polymers are neutral, $`R_g`$ is the only length scale appearing in this colloidal description, whereas if they carry charge, additional scales set, e.g., by the concentration of the solution, the counterions and/or the salt ions also come into play. We examine some characteristic cases below.

Polymer chains. Two possible choices for the effective coordinates have been investigated thus far. One possibility is to consider the effective interaction $`v_{\mathrm{com}}(r)`$ between the centres of mass of the linear chains, when these are kept fixed at a distance $`r`$ from one another. This was first done in the pioneering work of Flory and Krigbaum , who found a Gaussian interaction between the centers of mass. Although the functional form of the Flory-Krigbaum potential is correct, the dependence of the prefactor on the degree of polymerisation $`N`$ is not: whereas the Flory-Krigbaum mean-field approach predicts a $`N^{1/5}`$-dependence the prefactor turns out to be $`N`$-independent for sufficiently large $`N`$, as can be easily seen from examining the polymer second virial coefficient. Standard scaling arguments show that the radius of gyration $`R_g`$ is the only relevant length scale for the dilute and semi-dilute regimes of polymers in a good solvent . This immediately implies that the second-virial coefficient should scale as $`B_2R_g^3`$, which, in turn, implies that effective interaction must have an amplitude that is independent of $`R_g`$ , at least in the scaling limit.

A number of simulational  as well as theoretical approaches  involving two self-avoiding chains, reached the conclusion that the aforementioned interaction has a Gaussian form. The lack of divergence of this effective interaction at zero separation should not be surprising. Indeed, the centres of mass of two polymer chains can coincide without any of the monomers violating the excluded volume conditions. In addition, it can be seen that the effective “particles” one chooses for the coarse-grained description of the system do not need to be real particles of the physical system.

Recently, Louis et al.  have independently carried out state-of-the art simulations involving not just two but $`N_c`$ chains and varying the number of chains to cover a very broad range of concentrations, ranging from dilute solutions up to nine times the overlap concentration. They confirmed that the effective potential has a Gaussian-like form which at zero density can be well approximated by

$$v_{\mathrm{com}}(r)=\epsilon \mathrm{exp}[(r/\sigma )^2],$$
(2)
where $`\epsilon =1.87k_BT`$ and $`\sigma =1.08R_g`$. For higher densities a superposition of three Gaussians provides a very accurate fit , but the basic shape does not deviate much from the low-density Gaussian form. The rather weak density dependence can be shown to arise from the density independent many-body forces . Moreover, the same authors have shown that employing this effective interaction leads to a very accurate description of the thermodynamics (equation of state) of polymer solutions for a wide range of concentrations, thus confirming the validity of the idea that polymer chains can be viewed as soft colloids .

An alternative is to consider the end-monomers or the central monomers of the chains as effective coordinates . General scaling arguments establish that in this case the effective interaction diverges logarithmically with the monomer-monomer separation $`r`$ . When the central monomers are chosen, linear chains are equivalent to star polymers with $`f=2`$ arms. Motivated by this analogy, Jusufi et al.  derived the effective interaction $`v_{\mathrm{cm}}(r)`$ by combining monomer-resolved, off-lattice simulations with theoretical arguments. The sought-for interaction features in this case a logarithmic divergence for small separations, in full agreement with the scaling arguments, and crosses over to a Gaussian decay for larger ones:

$`v_{\mathrm{cm}}(r)={\displaystyle \frac{5}{18}}k_BTf^{3/2}\{\begin{array}{cc}\mathrm{ln}\left(\frac{r}{\sigma _\mathrm{s}}\right)+\frac{1}{2\tau ^2\sigma _\mathrm{s}^2}\hfill & \text{for }r\sigma _\mathrm{s}\text{;}\hfill \\ \frac{1}{2\tau ^2\sigma _\mathrm{s}^2}\mathrm{exp}\left[\tau ^2(r^2\sigma _\mathrm{s}^2)\right]\hfill & \text{for }r>\sigma _\mathrm{s}\text{,}\hfill \end{array}`$ (3)
where $`\sigma _\mathrm{s}0.66R_g`$ and $`\tau (f)`$ is a free parameter of the order of $`1/R_g`$ and is obtained by fitting to computer simulation results. For $`f=2`$ the value $`\tau \sigma _\mathrm{s}=1.03`$ has been obtained , which, together with the potential in Eq. (3) above yields for the second virial coefficient of polymer solutions the value $`B_2/R_\mathrm{g}^3=5.59`$, in agreement with the estimate $`5.5<B_2/R_\mathrm{g}^3<5.9`$ from renormalisation group and simulations .

Dendrimers. By employing a simple, mean-field theory based on the measured monomer density profiles of fourth-generation dendrimers, a Gaussian function of the form (2) has been shown to accurately describe the effective interaction between the centres of mass of these dendrimers . The prefactor $`\epsilon `$ has in this case a higher value then for linear polymers, $`\epsilon 10k_BT`$. Small-angle neutron scattering (SANS) profiles from concentrated dendrimer solutions are reproduced very well theoretically, at least below the overlap concentration $`c_{}`$.

Star polymers. By chemically anchoring $`f`$ linear chains on a common core, star polymers with functionality $`f`$ are constructed. In the theoretical analysis of the conformations and the effective interactions of stars, the finite size of the core particle is ignored, an excellent approximation when the chains are long. The natural choice for the effective coordinates is now the position of the central particle, i.e., of the star centre. For small functionalities, $`f10`$, Jusufi et al.  have shown that a logarithmic-Gauss potential of Eq. (3) accurately describes the effective interaction. The decay parameter $`\tau `$ of the Gaussian is $`f`$-dependent, for details see Ref. . For larger functionalities, $`f10`$, the Daoud-Cotton  blob picture of the stars is valid and the star-star interaction potential $`v_{\mathrm{ss}}(r)`$ reads as :

$`v_{\mathrm{ss}}(r)={\displaystyle \frac{5}{18}}k_BTf^{3/2}\{\begin{array}{cc}\mathrm{ln}\left(\frac{r}{\sigma _\mathrm{s}}\right)+\frac{1}{1+\sqrt{f}/2}\hfill & \text{for }r\sigma _\mathrm{s}\text{;}\hfill \\ \frac{\sigma _\mathrm{s}/r}{1+\sqrt{f}/2}\mathrm{exp}\left[\frac{\sqrt{f}}{2\sigma _\mathrm{s}}(r\sigma _\mathrm{s})\right]\hfill & \text{for }r>\sigma _\mathrm{s}\text{,}\hfill \end{array}`$ (4)
with the “corona diameter” $`\sigma _\mathrm{s}0.66R_g`$. Both star-star potentials, the one valid for $`f10`$, Eq. (3) and the one valid for $`f10`$, Eq. (4), show an ultrasoft logarithmic divergence as $`r0`$. The strength of the divergence is controlled by the functionality $`f`$, so that at the formal limit $`f\mathrm{}`$ the interaction (4) tends to the HS-potential.

Polyelectrolyte stars. If the polymer chains of a star polymer contain ionisable groups, the latter dissociate upon solution in a polar (aqueous) solvent, leaving behind charged monomers and resulting in a solution consisting of charged star polymers and counterions. The resulting macromolecules are called polyelectrolyte (PE) stars. In PE stars the chains are stretched due to the Coulomb repulsion of the charged monomers. The degree of stretching and condensation of counterions on the rods depends on the amount of charge and on the Bjerrum length. For moderate to high charging fractions, the effective interaction between the centres of the PE-stars have been analysed recently by means of computer simulations and theory . This interaction, $`v_{\mathrm{pes}}(r)`$, is dominated by the entropic contribution of the counterions that remain trapped within the star corona. For a broad range of functionalities and charge fractions, it can be accurately described by the fit:

$`{\displaystyle \frac{v_{\mathrm{pes}}(r)}{k_BT}}={\displaystyle \frac{\stackrel{~}{C}fN_c}{1\zeta }}\{\begin{array}{cc}\left[1\left(\frac{r}{\sigma }\right)^{1\zeta }\right]+\frac{2}{5}\left[\left(\frac{r}{\sigma }\right)^{2\zeta }1\right]+\frac{3(1\zeta )}{5(1+\kappa \sigma )}\hfill & \text{for }r\sigma \text{;}\hfill \\ \frac{3(1\zeta )}{5(1+\kappa \sigma )}\left(\frac{\sigma }{r}\right)\mathrm{exp}[\kappa (r\sigma )]\hfill & \text{for }r\sigma \text{,}\hfill \end{array}`$ (5)
where $`N_c`$ is the number of counterions, $`\sigma `$ the corona diameter of the PE-star, and $`\kappa `$ the inverse Debye screening length due to the free counterions. Finally, $`\stackrel{~}{C}`$ and $`\zeta `$ are fit parameters, where $`0<\zeta <1`$. The last condition ensures that the potential of Eq. (5) above tends to a finite value as $`r0`$. Hence, once more we are dealing with an ultrasoft interaction that varies slowly as the particle centres approach one another.

Ultrasoft interactions therefore describe quite a number of different systems. Their common characteristic is that the constituent particles are polymers of various architectures that dominate the spatial extent of the aggregates. In other words, one expects similar ultrasoft interactions to show up also when one deals with core-shell particles, consisting of a solid core and a polymeric shell, whenever the thickness of the latter greatly exceeds the radius of the former. In addition, the ultrasoft interactions can be tuned by controlling the number of arms, the charge, the length of the chains, the generation number (in the case of dendrimers) etc. Hence, it is useful to explore the general characteristics of this family of potentials and the ramifications they have on the structural and thermodynamic properties of the fluid and crystal phases of such systems.

## 3 Mean-field fluids

Motivated by the fact that effective interactions between polymeric colloids can be bounded (i.e., finite at all separations $`r`$), we examine here in general the properties of systems characterised by pair potentials of the form

$$v(r)=\epsilon \varphi (r/\sigma ),$$
(6)
with $`\varphi (x)<\mathrm{}`$ for all $`x`$. In Eq. (6) above, $`\epsilon `$ is an energy scale and $`\sigma `$ a length scale. Moreover, $`v(r)`$ is non-attractive, i.e., $`\mathrm{d}\varphi (x)/\mathrm{d}x0`$ everywhere. We introduce dimensionless measures of temperature and density as

$`t`$ $`=`$ $`{\displaystyle \frac{k_BT}{\epsilon }}=\left(\beta \epsilon \right)^1;`$ (7)
$`\eta `$ $`=`$ $`{\displaystyle \frac{\pi }{6}}\rho \sigma ^3={\displaystyle \frac{\pi }{6}}\overline{\rho },`$ (8)
where $`k_B`$ is Boltzmann’s constant and $`\rho =N/V`$ is the density of a system of $`N`$ particles in the volume $`V`$. We will refer to $`\eta `$ as the ‘packing fraction’ of the system.

The key idea for examining the high temperature and/or high density limit of such model systems in three and higher dimensions is the following. We consider in general a spatially modulated density profile $`\rho (𝐫)`$ which does not vary too rapidly on the scale $`\sigma `$ set by the interaction. At high densities, $`\rho \sigma ^31`$, the average interparticle distance $`a\rho ^{1/3}`$ becomes vanishingly small, and it holds $`a\sigma `$, i.e., the potential is extremely long-range. Every particle is simultaneously interacting with an enormous number of neighboring molecules and in the absence of short-range excluded volume interactions the excess free energy of the system  can be cast in the mean-field approximation (MFA) to be equal to the internal energy of the system :

$$F_{\mathrm{ex}}[\rho ]\frac{1}{2}\mathrm{d}^3r\mathrm{d}^3r^{}v(|𝐫𝐫^{}|)\rho (𝐫)\rho (𝐫^{}),$$
(9)
with the approximation becoming more accurate with increasing density. Then, Eq. (9) immediately implies that in this limit the direct correlation function $`c(r)`$ of the system, defined as

$$c(|𝐫𝐫^{}|;\rho )=\underset{\rho (𝐫)\rho }{lim}\frac{\delta ^2\beta F_{\mathrm{ex}}[\rho (𝐫)]}{\delta \rho (𝐫)\delta \rho (𝐫^{})},$$
(10)
becomes independent of the density and is simply proportional to the interaction, namely

$$c(r)=\beta v(r).$$
(11)
Using the last equation, together with the Ornstein-Zernike relation , we readily obtain an analytic expression for the structure factor $`S(k)`$ of the system as

$$S(k)=\frac{1}{1+\overline{\rho }t^1\stackrel{~}{\varphi }(k\sigma )},$$
(12)
where $`\stackrel{~}{\varphi }(q)=\mathrm{d}^3x\mathrm{exp}(\mathrm{i}𝐪𝐱)\varphi (x)`$ is the Fourier transform of $`\varphi (x)`$.

Bounded and positive-definite interactions have been studied in the late 1970s by Grewe and Klein . The authors considered a Kac potential of the form:

$$v(r)=\gamma ^d\psi (\gamma r),$$
(13)
where $`d`$ is the dimension of the space and $`\gamma 0`$ is a parameter controlling the range and strength of the potential. Moreover, $`\psi (x)`$ is a nonnegative, bounded and integrable function. Grewe and Klein showed rigorously that at the limit $`\gamma 0`$, the direct correlation function of a system interacting by means of the potential (13) is given by Eq. (11) above. The connection with the case we are discussing here is straightforward: as there are no hard cores in the system or a lattice constant to impose a length scale, the only relevant length is set by the density and is equal to $`\rho ^{1/3}`$ in our model and by the parameter $`\gamma ^1`$ in model (13). In this respect, the limit $`\gamma 0`$ in the Kac model of Grewe and Klein is equivalent to the limit $`\rho \mathrm{}`$ considered here.

Although the limit of Grewe and Klein corresponds to $`t\mathrm{}`$ and $`\overline{\rho }\mathrm{}`$, the relation (11) has been shown to be an excellent approximation at arbitrarily low temperatures for high enough densities  and for temperatures $`t1`$ practically at all densities . Hence, the mean-field approximation is valid in a vast range of the thermodynamic space of such systems, which has led to their characterisation as mean-field fluids (MFF) . Associated with the structural relation (11) are scaling relations of thermodynamic quantities, arising from the compressibility sum rule :

$$f^{\prime \prime }(\rho )=\stackrel{~}{c}(k=0;\rho )=4\pi _0^{\mathrm{}}r^2c(r;\rho )dr,$$
(14)
where $`f(\rho )=\beta F_{\mathrm{ex}}(\rho )/V`$, and the primes denote the second derivative. From Eqs. (11) and (14) it then follows

$$f(\rho )=\frac{\beta \stackrel{~}{v}(0)}{2}\rho ^2,$$
(15)
with $`\stackrel{~}{v}(0)=4\pi _0^{\mathrm{}}r^2v(r)dr`$. This simple scaling is not at all equivalent to a second virial theory. In fact, simple virial expansions have a rather small radius of convergence for mean field fluids . It then follows that the excess pressure, chemical potential and compressibility satisfy the scaling relations $`P_{\mathrm{ex}}\rho ^2`$, $`\mu _{\mathrm{ex}}\rho `$ and $`\chi _{\mathrm{ex}}\rho ^2`$ .<sup>1</sup><sup>1</sup>1Note that for polymers in a good solvent $`P_{\mathrm{ex}}\rho ^{3\nu /(3\nu 1)}\rho ^{2.3}`$ in the semi-dilute regime. The density dependence of the pair-potentials is necessary to properly describe this correction to simple MFF behaviour . But this in turn implies that the extra factor $`0.3`$ in the scaling of the pressure arises from the many-body interactions, since these are what cause the density-dependence in the first place . All these stem from the validity of the strong structural-MFA relation (11) which guarantees the validity of the weaker thermodynamical MFA relation (15). In what follows, we will argue that many ultrasoft systems can still satisfy the thermodynamic relation (15) approximately, without satisfying Eq. (11). To distinguish between the two classes, we now qualify the term “mean-field fluids” and call the ones for which both the structural and thermodynamic MFA work well strong mean-field fluids. For such systems, for which the potentials are typically also bounded, the density functional of Eq. (9) has been extended to mixtures  allowing a straightforward and transparent analysis of the phase separation in the bulk, interfacial  and wetting properties  of such mixtures.

Let us now then turn our attention to interaction potentials such as those given in Eqs. (3) and (4). These diverge at the origin slowly enough, so that the three-dimensional integral

$$\mathrm{d}^3rv(r)=4\pi _0^{\mathrm{}}r^2v(r)dr=\stackrel{~}{v}(0),$$
(16)
is finite and equal (by definition) to the value of the Fourier transform of the potential $`\stackrel{~}{v}(k)`$ at $`k=0`$. It is now impossible to satisfy the strong mean-field condition of Eq. (11) everywhere. Indeed, the direct correlation function $`c(r)`$ has to remain finite at $`r=0`$, whereas the pair potential diverges. Hence, as shown in Fig. 1(a), there will always exist a region in the neighborhood of the origin in which Eq. (11) is violated. At the same time, it can be seen in this figure that the extent of this region shrinks with increasing density, hence the fluid becomes more ‘strong mean-field’-like as it gets denser. The smaller $`\stackrel{~}{v}(0)`$, the lower the density at which the MFA for $`S(k)`$, Eq. (12), becomes a reasonable approximation.

The discrepancies between $`c(r)`$ and $`\beta v(r)`$ become much less important when we turn our attention to the thermodynamics. To obtain the excess Helmholtz free energy, one needs only the integral of $`r^2c(r)`$, see Eq. (14). As demonstrated in Fig. 1(b), upon multiplication with the geometrical factor $`r^2`$, the deviations of $`c(r)`$ from $`\beta v(r)`$ become suppressed, so that we can write, to a very good approximation:

$$_0^{\mathrm{}}drr^2c(r;\rho )_0^{\mathrm{}}drr^2\beta v(r).$$
(17)
Eq. (17) together with Eq. (14) yield an approximate scaling of the excess free energy of the weak mean-field fluids with density that is identical to that of the strong mean-field fluids, Eq. (15). The accuracy of the approximation for the star polymer fluid with $`f=32`$ arms \[Eq. (4)\] is shown in Fig. 2. The line labeled as exact free energy there was obtained by solving the Rogers-Young closure  for the fluid at a wide density range and subsequently utilizing the compressibility sum rule \[Eq. (14)\] to obtain the excess free energy. Comparisons with simulations  have indeed demonstrated that this procedure delivers an essentially exact numerical result.

Clearly, the mean-field approximation improves with increasing density, as the number of particles effectively interacting with one another grows. The crossover density $`\rho _\times `$ above which the quadratic scaling of the free energy holds is $`f`$-dependent and grows with increasing $`f`$. Indeed, the functionality acts as a prefactor that controls the strength of the logarithmic divergence of the potential at the origin. Formally, the mean-field approximation also becomes better with growing spatial dimension $`d`$, as the geometrical prefactor $`r^{d1}`$ multiplying $`c(r)`$ and $`\beta v(r)`$ suppresses the small-$`r`$ discrepancies of the two more efficiently. We call systems for which the mean-field idea holds only for the thermodynamics weak mean-field fluids. If one naïvely applies the strong mean-field relation, Eq. (11), to weak mean-field fluids, one obtains results for the structure factor $`S(k)=[1\rho \stackrel{~}{c}(k)]^1`$ that can be seriously in error for finite $`k`$-values. Only at $`k=0`$ and at sufficiently high densities is it a reasonable approximation to set $`S(0)=[1+\rho \beta \stackrel{~}{v}(0)]^1`$.

## 4 Clustering and reentrant melting

In this section, we turn our attention to the phase behaviour of strong mean-field fluids. Two representatives of this class whose phase behaviour has been studied in detail are the Gaussian core model (GCM) of Eq. (2) and the ‘penetrable sphere model’ (PSM) characterised by the interaction potential $`v_{\mathrm{psm}}(r)=\epsilon \mathrm{\Theta }(\sigma r)`$, with the Heaviside step function $`\mathrm{\Theta }(x)`$. Clearly, the PSM reduces to the hard sphere model for $`t=0`$.

The GCM has been the subject of extensive investigations by Stillinger et al. in the late 1970s . The $`t=0`$ phase diagram of the model was calculated, showing the existence of two stable crystal structures, fcc for low densities and bcc for high ones. In addition, a host of mathematical relations for the GCM has been established and on the basis of free energy estimates is has been postulated that the system displays reentrant melting behaviour at low temperatures. On the basis of simulation studies at selected thermodynamic points, a rough phase diagram of the GCM has been drawn . A detailed study of the structural and phase behaviour of the GCM was carried out recently by Lang et al. . There, it was indeed shown that for temperatures $`t>0.01`$ the system remains fluid at all densities, whereas for $`t0.01`$ reentrant melting is observed: increasing the density, the system first undergoes a fluid $``$ fcc transition, followed by a structural fcc $``$ bcc transition and at higher densities the bcc solid remelts, i.e., a bcc $``$ fluid transition takes place. The width of the solid-phase region grows with decreasing temperature. A structural signature of this unusual phase diagram in the fluid phase above $`t=0.01`$ is an anomaly in the behaviour of the liquid structure factor $`S(k)`$. The height of its main peak first grows with increasing density and after achieving a maximum, it decreases again, reflecting the stability of the fluid beyond the reentrant melting. The Hansen-Verlet freezing criterion  was shown to be satisfied at both the freezing and the reentrant melting lines . Moreover, it was found that at high densities not only the mean-field relation, Eq. (11) is satisfied to excellent accuracy but also that the hypernetted chain closure (HNC) becomes quasi-exact . The system becomes ‘quasi-ideal’ at those densities, meaning that the radial distribution function has the limiting behaviour $`g(r)1`$. This, in conjunction with the mean-field property $`c(r)=\beta v(r)`$ and the exact relation $`g(r)=\mathrm{exp}[\beta v(r)c(r)+g(r)1+B(r)]`$, forces the bridge function to obey the limit $`B(r)0`$, hence rendering the HNC exact.

The PSM was first studied in detail by by means of integral equation theories, computer simulations and cell-model calculations at small temperatures, $`t0.3`$ . In contrast to the GCM, no reentrant melting was found. Instead, the freezing line of the system appeared to persist at all temperatures, and cascades of clustering transitions in the solid were found, in which solids with multiple site occupancies are stable with increasing temperature and density. These findings were independently confirmed in a density-functional study of the low-temperature phase behaviour of the PSM . Sophisticated integral-equation approaches at arbitrarily high temperatures revealed a loss of the solution along the ‘diagonal’ $`t=\eta `$ of the phase diagram , again a feature pointing to an instability of the liquid by increasing density at arbitrarily high temperatures. Thus, the PSM and the GCM show completely different phase behaviours, although they are both bounded and non-attractive potentials. There is a cascade of clustering transitions for the former, enabling freezing at all temperatures, and a reentrant melting for the latter, associated with the inability to stabilise crystals above a certain critical temperature.<sup>2</sup><sup>2</sup>2The use of the terminology ‘critical temperature’ here refers simply to the fact that above the said temperature freezing is impossible and should not be confused with its standard meaning in the realm of critical phenomena. There are no diverging thermodynamic quantities here, no universality and all free energies remain analytical functions in the neighborhood of the critical temperature.

The key in understanding these two very different types of behaviour lies in the strong mean-field character of these fluids and the associated expression for the fluid structure factor, Eq. (12). If the Fourier transform of the pair interaction $`\stackrel{~}{\varphi }(k)`$ has oscillatory behaviour (i.e., if $`\stackrel{~}{\varphi }(k)`$ becomes negative for some $`k`$-values), then at the wavenumber $`k_{}\sigma `$ where $`\stackrel{~}{\varphi }(k\sigma )`$ attains its most negative value, $`|\stackrel{~}{\varphi }(k_{}\sigma )|`$, the liquid structure factor $`S(k_{})`$ will display a maximum. For any given temperature $`t`$, there exists then a density $`\overline{\rho }_\mathrm{s}(t)`$ such that $`\overline{\rho }_\mathrm{s}|\stackrel{~}{\varphi }(k_{}\sigma )|=t`$, causing a divergence of the fluid structure factor at $`k_{}`$ and marking a ‘spinodal line’ at this finite wavenumber. Thus, the fluid cannot be stable at all densities. As a matter of fact, freezing will take place before the spinodal line is reached. The PSM clearly belongs to this category, since the abrupt jump of the pair interaction $`v_{\mathrm{psm}}(r)`$ at $`r=\sigma `$ causes long-range oscillations of the potential in Fourier space. By employing the Hansen-Verlet criterion, $`S(k_{})3`$, it has been found  that in the PSM freezing takes place at the ‘diagonal’ $`t=\eta `$ on which the integral equation approach of Fernaud et al.  breaks down. If, on the other hand, the Fourier transform of the potential, $`\stackrel{~}{\varphi }(k)`$ is a positive-definite, monotonically decreasing function of $`k`$, Eq. (12) assures that at sufficiently high temperatures, where the mean-field approximation is valid at all densities , $`S(k)`$ is a monotonic function of $`k`$ approaching rapidly the value $`S(k)=1`$ with increasing $`k`$ and being deprived of any peaks. The lack of peaks in the structure factor implies the lack of any tendency within the liquid towards spontaneous formation of spatially modulated patterns. Thus, the fluid remains stable for all densities at sufficiently high temperatures. This, combined with the observation that at low temperatures and densities bounded potentials all become hard-sphere like and hence they must cause a freezing transition there, leads to a reentrant-melting scenario for such systems. Clearly, the GCM belongs to this category. Representative results for a particular family of strong mean-field fluids and schematic phase diagrams can be found in Ref. .

## 5 Exotic crystal phases

We now focus on weak-mean field fluids, for which no simple criterion for their freezing behaviour can be established, since Eqs. (11) and (12) are not satisfied any more. The star-polymer fluid characterised by the pair potential of Eq. (4) is a case in point.<sup>3</sup><sup>3</sup>3We consider here the case $`f>10`$, for which indeed the interaction of Eq. (4) holds, and not the low-functionality case for which the interaction of Eq. (3) is valid. The reason is that at low functionalities the stars do not freeze at any density and hence they are not an appropriate system for considering thermodynamically stable crystals. The physical system of star polymers provides an excellent testbed for the investigation of the thermodynamic stability of more complicated crystals than usual fcc- and bcc-lattice arrangements.

The fcc-lattice is the one favoured by hard interactions, since it has the property of maximising the available volume and hence the entropy of the particles for a given particle density . On the other hand, the presence of ‘soft tails’ in the potential has the effect of favouring the more open bcc-lattice, as was convincingly demonstrated for the case of the screened Coulomb potential (Yukawa interaction) arising in charge-stabilised colloidal suspensions . These two common lattices were considered for a long time to be the only ‘candidates’ in a search for stable crystals for given interatomic potentials. However, in modern colloidal science, new possibilities open up. It is technically possible to manufacture micelle-like particles featuring a hard core and a soft, fluffy corona of grafted or adsorbed polymer chains, with the thickness $`L`$ of the latter being much larger than the radius $`R_c`$ of the former. Star polymers correspond to the case $`R_cL`$; the theoretical arguments leading to the effective potential of Eq. (4) are in fact based on the assumption $`R_c0`$. Under these physical circumstances, the effective interaction between the micelles is dominated by the ultrasoft repulsion between the overlapping, flexible coronas and not by the excluded volume interactions of the hard cores. Thereby, the requirement of maximising of the volume available to each particle does not play the decisive role any more.

These properties manifest themselves in the phase diagram of star-polymer solutions. Due to the irrelevance of the temperature for this entropic interactions, the phase diagram was drawn in the ($`f,\eta `$)-plane by Watzlawek et al. , where $`\eta =\pi \rho \sigma ^3/6`$ and $`\rho =N/V`$ is the number density of $`N`$ stars in the volume $`V`$. The phase diagram is shown in Fig. 3. The fluid phase remains stable at all concentrations for $`f<f_c=34`$, a result that confirms and makes precise early scaling-argument predictions of Witten et al. . For $`f>f_c`$ and at packing fractions $`0.15\eta 0.70`$, the usual fcc- and bcc-crystals are seen to be stable, the former for larger and the latter for smaller functionalities. This is consistent with the fact that the effective potential of Eq. (4) has a Yukawa decay length scaling as $`1/\sqrt{f}`$, hence large $`f`$ is analogous to the strongly screened charge-stabilised colloidal suspensions. However, for $`\eta 0.70`$, unusual crystal structures appear. First, in the domain $`0.70\eta 1.10`$, a body-center-orthogonal (bco) crystal is thermodynamically stable. The bco-lattice is characterised by a body-centered, orthogonal conventional unit cell with three unequal sides and reduces to the bcc-lattice for ratios $`1:1:1`$ between the sides and to the fcc for ratios $`1:1/\sqrt{2}:1/\sqrt{2}`$ . The bco-lattices appearing in this region of the phase diagram feature strongly anisotropic unit cells with typical size ratios $`1:0.6:0.3`$. Thus, these are crystals with coordination number 2. For packing fractions $`1.10\eta 1.50`$, the diamond lattice with coordination number 4 turns out to be stable. Thus, we see that very open structures with their characteristically low coordination numbers are stabilised by the ultrasoft star-star potential. This feature has been attributed to the very slow divergence of the interaction as $`r0`$, combined with its crossover to a Yukawa-form for $`r>\sigma `$ . Indeed, in such circumstances it may be energetically preferable for the system to have a small number of nearest neighbours at a small distance from any given lattice point than a large number of neighbours at a greater distance, as is the case for the optimally-packed fcc-lattice.

In the fluid there exist clear structural signatures both for the topology of the phase diagram of Fig. 3 and for the variety of the crystal phases featured there. The reentrant transition from a fluid to a bcc-lattice and then again to a fluid, occurring for $`34f48`$, is manifested in fluid structure factors $`S(q)`$ that show a main peak that first grows with increasing density and then drops again , as in the case of the Gaussian core model mentioned above. Moreover, once more the Hansen-Verlet freezing criterion  was found to be satisfied on both sides of the freezing and reentrant melting line. The radial distribution function $`g(r)`$ of the fluid at various densities, on the other hand, carries the signature of a local coordination that resembles that of the thermodynamically neighbouring solids. To demonstrate this, we show in Fig. 4 the function $`g(r)`$ for star polymer fluids at $`f=32`$, which are thermodynamically stable, at packing fractions $`\eta =0.80`$ and $`\eta =1.20`$. Comparison with Fig. 3 shows that the former corresponds to a state at the vicinity of the bco-phase and the latter to one at the vicinity of the diamond phase. Consider now the average coordination number $`z`$ in the fluid phase, defined as

$$z=4\pi \rho _0^{r_{\mathrm{min}}}r^2g(r)dr,$$
(18)
where $`r_{\mathrm{min}}`$ is the position for which $`g(r)`$ has its first minimum and is indicated by the arrows in Fig. 4. For the two packing fractions shown we obtain $`z=1.95`$ at $`\eta =0.80`$ and $`z=4.03`$ at $`\eta =1.20`$. The first is very close to the coordination number $`z_{\mathrm{bco}}=2`$ of the neighboring bco-lattice and the latter to $`z_{\mathrm{diam}}=4`$ of the diamond lattice. The fluid distribution functions contain local correlations that point to the ordering of the incipient crystal phases.

The $`g(r)`$’s of the fluid above the overlap density, $`\eta 1.0`$, show in addition anomalous behaviour featuring two distinct length scales, as is clear from Fig. 4(b). As analysed in detail in Ref. , two characteristics of the interaction are responsible for this behaviour: on the one hand, the existence of the crossover of the interaction of Eq. (4) at $`r=\sigma `$ from a logarithmic to an exponentially decaying form. And on the other, the ultrasoftness of the logarithmic potential, allowing the existence of fluids at arbitrarily large densities (for $`f<f_c`$), a feature unknown for the usual interactions appearing in liquid-state physics and which are all ‘perturbations’ of the Hard-Sphere potential (e.g., Lennard-Jones, inverse-powers etc.) Thus, ultrasoft potentials carry unique structural signatures which should in principle be visible in scattering experiments from soft, polymeric fluids.

A great deal of insight into the general physical mechanisms driving the stability of open structures in soft systems was gained through the recent work of Ziherl and Kamien . They considered in full generality systems with particles composed of a hard core and long, deformable coronas and argued as follows. At any given density above the overlap concentration, the coronas are forced to overlap and compress, which gives rise to an entropic free energy cost. The volume available to the coronas is fixed and equal to the difference of the total volume minus that occupied by the hard cores. Denoting by $`d`$ the thickness of the coronal layer and $`A`$ the total area of an imaginary membrane separating the compressed coronas, it then turns out that the product $`Ad`$ is constant. As the free energy cost for the compression of the chains increases with decreasing thickness $`d`$, it turns out that favourable phases are those for which the interfacial area $`A`$ is minimal. Thereby, the problem reduces to that of determining the ordered arrangement of point particles that generates Wigner-Seitz (WS) cells having the smallest possible area for a given density. It is then conceivable that the fcc-lattice will be unfavoured, since its WS cell has a larger area than that of the bcc, for instance. In this way, Ziherl and Kamien established a beautiful connection of this problem with Lord Kelvin’s celebrated question of determining the area-minimising partition of space for an arrangement of soap bubbles of equal volume .

Following these arguments, it then turns out that there exists yet another candidate phase that has an area even smaller than the bcc-lattice , namely the A15-lattice . Self-assembled micelles of dendritic molecules with a particular architecture have been experimentally seen to crystallise into this phase . The conventional unit cell of the A15-lattice is shown in Fig. 5. It can be thought of as the cell of a bcc-lattice (dark points) decorated with ‘dimers’ (light points) running along the middle of the faces. The dimers are oriented parallel on opposite phases, forming thus columns through space. The orientation of the dimers lying on intersecting faces is perpendicular to one another, so that one third of all dimers lies along each of the three Cartesian directions in space. The dimer length is $`a/2`$, where $`a`$ is the edge length of the cube, and it is placed symmetrically along the face, i.e., the distance of any monomer to its nearest edge is $`a/4`$. The A15-lattice is not a Bravais lattice; it can be constructed as a simple cubic (sc) arrangement with an eight-member basis, thus it contains 8 sites per conventional cell. Its WS-cell is a Goldberg decatetrahedron<sup>4</sup><sup>4</sup>4We prefer the term decatetrahedron for a polyhedron with 14 faces instead of the term tetrakaidecahedron, often employed in the literature. The latter, inspired from ancient Greek, literally means ‘four-and-ten-faced polyhedron’, whereas the former, consistent with modern Greek, has the much more logical translation ‘fourteen-faced polyhedron’. One of us (CNL) believes that modern Greek language deserves a fair chance against its classical predecessor, which has influenced scientific terminology for quite some time. consisting of two hexagonal and 12 pentagonal faces .

Semi-quantitative calculations on the stability of the A15-lattice were carried out by Ziherl and Kamien , using a model ‘square-shoulder’ potential within a simplified cell model. Narrow regions in thermodynamic phase space were found, in which the A15-lattice was stable but this finding is uncertain in view of the approximations involved and the limited extent of this region. To investigate this question in more detail, we have employed extensive lattice-sum calculations for the star-polymer system, using the pair potential of Eq. (4) and extending both the set of candidate lattices and the region of densities we looked at. We compared between the sc, diamond, bco (which includes the fcc- and bcc-lattices as special cases) and A15-lattices in the regions $`0\eta 2.50`$ and $`32f256`$, at selected arm numbers $`f`$ to be shown below. The lattice sums were performed by fixing the particles at the prescribed lattice positions and keeping them frozen there, i.e., no thermal fluctuations (harmonic corrections) were taken into account. This approach reproduces very well the solid part of the phase diagram of Fig. 3: the free energy of the star-polymer crystals turns out to be dominated by the lattice-sum term, the corrections to it from the particle oscillations as well as the entropic contribution from the same playing only a minor role. In this way, we are of course unable to compare the solid free energies with those of the fluid, therefore no prediction about melting can be made. However, for large enough $`f`$, the interaction is steep enough, so that the system will be definitely be in a crystalline phase for which the lattice sums provide an accurate prediction of the most stable structure among the candidates.

For the bco-phases we minimised the lattice sums with respect to the two size ratios, $`r_1=b/a`$ and $`r_2=c/a`$ between the edge lengths, $`a`$, $`b`$ and $`c`$ at any given density. Without loss of generality, we assume in what follows that $`a`$ is the longest of the three edges, thus $`0<r_1,r_21`$. The minimised bco-energy was then compared with the energies of all other lattices; the one with the smallest lattice energy per particle wins.

The ‘zero-temperature’ phase diagram<sup>5</sup><sup>5</sup>5‘Zero-temperature’ is just a convenient way to refer to the assumption of frozen particles at the lattice sites and does not refer to the real temperature $`T`$ of the system. The latter is an irrelevant thermodynamic variable because the entropic interaction of Eq. (4) is proportional to $`k_BT`$ and the thermal energy is the only energy scale of the problem. obtained this way is shown in Fig. 6. First, we note that the phases being stable up to $`\eta 1.50`$ are precisely those also seen in the finite-temperature phase diagram of Fig. 3 and also that the phase boundaries based on the lattice sums agree very well with the ones at finite temperatures. The A15-phase does not alter the hitherto explored part of the phase diagram. However, for $`\eta 1.50`$, a host of new phases and of transitions between those show up. The A15-phase is stable at the high-density part of the phase diagram, confirming thus explicitly the prediction of Refs.  and  that this phase is a suitable candidate at the high concentrations of ultrasoft particles. Nested between the stability domain of the A15-lattice and the diamond lattice, four new bco-phases and (iso)structural transitions between those show up, having the following characteristics.

The phase denoted bco<sub>2</sub> has size ratios $`r_1=1`$ and $`r_20.55`$, the former being constant at all densities and functionalities and the latter showing very weak variation, see also Fig. 7. Hence, the unit cell of the bco<sub>2</sub>-phase is anisotropic only in one Cartesian direction and has a wide basis and height that is smaller than the base edge length. Accordingly, the coordination number of this phase is 2. The bco<sub>3</sub>-phase, dominating at high functionalities, has size ratios $`r_1=r_20.3`$ which also show very little variation with $`\eta `$ and $`f`$, see Fig. 7. Similarly to the bco<sub>2</sub>-phase, therefore, the cell of the bco<sub>3</sub>-phase also has anisotropy in only one Cartesian direction. Contrary to it, however, the height of the cell is now larger than the edge length of the base and therefore the number of nearest neighbours is 4. The phase transition bco<sub>2</sub> $``$ bco<sub>3</sub> is first order: as can be seen in Figs. 7(c) and (d), the size ratios jump abruptly from the values ($`1,0.55`$) of bco<sub>2</sub> to the values ($`0.3,0.3`$) of bco<sub>3</sub>. These are the only stable bco-phases in this part of the phase diagram as long as $`f56`$.

For functionalities $`32f56`$, two more bco-phases show up. First, there is a narrow region occupied by the bco<sub>4</sub>-phase, which has all three edge lengths of its unit cell different and is thus similar to the bco<sub>1</sub>-phase discovered before. As seen in Figs. 7(a), (b) and (c), these ratios evolve from the values $`r_1=1`$ and $`r_20.55`$ of the bco<sub>2</sub>-phase towards the values $`r_1=1/\sqrt{2}`$ and $`r_2=1/\sqrt{2}`$ of the fcc-phase.

The order of the transitions between those phases has been investigated numerically. Within the limits of accuracy of our numerical code, we have found that for $`f=32`$ and $`f=40`$ the transition bco<sub>4</sub> $``$ fcc is first order, i.e., the size ratios of the bco<sub>4</sub>-phase jump with density to the ratios $`r_1=r_2=1/\sqrt{2}`$ of the fcc-phase abruptly. This is seen in Figs. 7(a) and Fig. 7(b). For $`f=48`$, the transition is second-order with the bco<sub>4</sub>-size ratios evolving to the fcc-ones smoothly, see Fig. 7(c). The nature of the transition bco<sub>2</sub> $``$ bco<sub>4</sub> is also mixed. For $`f=32`$, \[Fig. 7(a)\], a finite jump of the values of the ratios was found, although the step size in changing the packing fraction was made as small as $`5\times 10^5`$. Thus, we characterise this transition at $`f=32`$ as first-order. For $`f=40`$, \[Fig. 7(b)\] this transition appears to be second-order. Hence, we conclude that there must be a line of second-order transitions terminating at a tricritical point between $`f=40`$ and $`f=32`$, (for the bco<sub>2</sub> $``$ bco<sub>4</sub> transition) and similarly a tricritical point between $`f=48`$ and $`f=40`$ (for the bco<sub>4</sub> $``$ fcc-transition) to be succeeded by lines of first-order transitions.

All four bco-phases and the three transitions among them, bco<sub>2</sub> $``$ bco<sub>3</sub> $``$ bco<sub>4</sub> $``$ fcc, appear only in a very narrow $`f`$-range around $`f=48`$, see Fig. 7(c). The line of second-order transitions bco<sub>2</sub> $``$ bco<sub>4</sub> meets the lines of first-order transitions bco<sub>2</sub> $``$ bco<sub>3</sub> and bco<sub>3</sub> $``$ bco<sub>4</sub> at a critical endpoint located at $`f`$ slightly less than 48. Similarly, the line of second-order transitions bco<sub>4</sub> $``$ fcc also terminates at a critical endpoint, meeting the first-order lines bco<sub>3</sub> $``$ A15 and fcc $``$ A15. It is an intriguing phenomenon that such a richness in the stable crystal structures and in the nature of the transitions among them occurs as the result of a simple, spherically symmetric interaction and this points to the many surprises of ultrasoft potentials and their tendency to produce open, exotic structures. At the same time, it must be pointed out that we do not expect the solid phases occurring for $`f<48`$ to survive the competition with a fluid. Indeed, as can be seen from Fig. 3, the critical functionality $`f_c`$ below which no solids are stable increases with increasing density. The bcc-phase is extinguished by the fluid for $`f<34`$, whereas the bco- and diamond phases for $`f<44`$. It is therefore anticipated that the bco<sub>2</sub>-, bco<sub>4</sub>\- and fcc-phases seen in Fig. 6 for $`f<48`$ will be wiped out by the fluid there. However, by arbitrarily increasing $`f`$ one can always reach a domain where the fluid will be beaten by the crystal and hence the bco<sub>2</sub> $``$ bco<sub>3</sub> $``$ A15-transitions will be there also after a finite-temperature calculation. In particular, the A15-phase, on whose stability it has been speculated for some time, has now been proved to be indeed the most stable phase at sufficiently high densities among the candidates considered.

The logarithmic-Yukawa potential employed in our study has been shown to describe well the effective interactions between star polymers in a good solvent for a wide range of concentrations. However, at the region of stability of the A15-crystal, the pair potential description is not expected to be particularly accurate. Many-body contributions are expected to become important there . Moreover, the Yukawa tail of the potential, describing the interactions of the outermost Daoud-Cotton blobs, should be absent since the compressed coronas there are deprived of the outermost blob structure of the isolated stars. Thus, in this respect, the logarithmic-Yukawa potential has to be looked upon rather as a toy model. Nevertheless, the physical characteristic driving the transitions discovered above is the ultraslow divergence of the logarithm which, in the neighborhood of the average particle distance can be locally expanded as a ramp-like potential and the finer details of the interaction should become irrelevant.

## 6 Summary and conclusions

We have shown that “ultrasoft interactions” arise naturally from coarse-graining procedures for a broad range of soft matter systems. Besides greatly simplifying the statistical mechanics of these complex systems – once the interactions are derived, all the well known tools of liquid state theory can be applied to calculate correlations and phase behaviour – they also lead to new phenomenology. Signatures of these ultra soft interactions include anomalous fluid correlations, reentrant melting as well as the stabilisation of exotic, open crystal structures. In contrast to their atomic counterparts, soft matter systems can therefore stabilise such crystals without the presence of angle-dependent, anisotropic potentials: radially symmetric, ultrasoft interactions are quite sufficient. Thus, a new mean field fluid paradigm can be established, which goes beyond the usual prototype for classical fluids, the hard-sphere model. The latter, being always dominated by packing effects, tends to favour close packed structures. Exotic crystals with unusual ordering have been observed in hard-sphere-like suspensions  but that case refers to binary mixtures whose phase behaviour is indeed much richer than that of their one-component counterparts.

It is with great pleasure that we dedicate this paper to Professor Peter Pusey on the occasion of his 60th birthday. We thank Martin Watzlawek and Primoz Ziherl for helpful discussions. AAL thanks the Isaac Newton Trust, Cambridge, for financial support.

## References


text:
# Ultrahigh energy cosmic rays from dark matter annihilation

## 1 INTRODUCTION

Explaining the origin of cosmic rays at energies above the Greisen–Zatsepin–Kuzmin bound $`E_{GZK}4\times 10^{19}`$eV is an interesting and puzzling scientific challenge for several reasons:

AGASA observes a flux of order of one cosmic ray event per $`100\mathrm{km}^2`$ per year at energies above $`10^{20}`$eV , thus confirming the flux seen by previous smaller ground based arrays on the northern hemisphere111But see  for a recent discussion of the pertinent problem to reconstruct primary cosmic ray energies from air showers. A new analysis presented there lowered the estimate on the energy of the highest energy event observed with the Haverah Park Array from above $`10^{20}`$eV to $`8.3\times 10^{19}`$eV and also indicates a smaller flux at the high energy end of the spectrum. Of course, the very existence of UHECRs is not questioned by these findings., see e.g. .

On the other hand, energy loss of electromagnetically interacting particles due to scattering with the CMB background indicates that these ultrahigh energy cosmic rays (UHECRs) should originate within a radius of order $`100`$Mpc , and sufficiently powerful astrophysical acceleration mechanisms are not easy to find.

This has motivated proposals that superheavy dark matter may be the origin of ultrahigh energy cosmic rays:

This includes in particular Hill’s proposal of monopolonium decay , the decay of superheavy dark matter relics proposed by Berezinsky et al.  and by Kuzmin and Rubakov , and most recently our proposal of collisional annihilation of superheavy dark matter . While decay scenarios for UHECRs should imply domination of the galactic center in the distribution of arrival directions, the proposal of collisional annihilation predicts that with sufficient statistics from the upcoming Auger observatory a clumpy structure should emerge in the galactic source distribution: Unitarity bounds on reaction cross sections imply that a smooth superheavy dark matter distribution in the galactic halo would not yield the observed flux through collisional annihilation, but the overdensity in dark matter clumps in the halo may well account for the observed flux .

The common virtues of these top-down scenarios are:
$``$ a natural explanation of the seemingly random distribution of arrival directions;
$``$ avoidance of collisional damping and the GZK bound due to origination in the galactic halo;
$``$ no need for an extremely powerful astrophysical acceleration mechanism.

Very high-energetic neutrinos both from superheavy dark matter decay and from annihilation have very recently been discussed in , and some recent overviews of top-down scenarios for UHE cosmic rays can be found in .

Superheavy dark matter particles cannot arise as a consequence of standard thermal freeze out, but it has been pointed out that dark particles might have never experienced local chemical equilibrium during the evolution of the universe, and that their mass may be in the range $`10^{12}`$ to $`10^{19}`$ GeV, much larger than the mass of thermal wimps . Since these wimps would be much more massive than thermal wimps, such superheavy dark particles have been called wimpzillas .

wimpzillas may be created during bubble collisions if inflation is completed through a first-order phase transition ; at the preheating stage after the end of inflation with masses easily up to the GUT scale  or even up to the Planck scale ; or during the reheating stage after inflation . wimpzillas may also be generated in the transition between an inflationary and a matter-dominated (or radiation-dominated) universe due to the “nonadiabatic” expansion of the background spacetime acting on the vacuum quantum fluctuations . The distinguishing feature of this mechanism is the capability of generating particles with mass of the order of the inflaton mass even when the particles only interact extremely weakly (or not at all) with other particles.

The mass $`m`$ of unstable dark matter which decays through gravitational couplings is constrained by the fact that the lifetime is only of order $`\tau m_{Planck}^2/m^3`$ , and therefore unstable heavy relics must decay through strongly suppressed channels which effectively correspond to sub-gravitational couplings. This provided one motivation for the collisional annihilation proposal in , since the need for 2-particle collisions implies a larger lifetime of annihilating wimpzillas as opposed to decaying superheavy particles.

## 2 THE FLUX FROM ANNIHILATION IN DARK MATTER CLUMPS

The spectral fluxes from decay or annihilation of wimpzillas of mass $`M_X`$ and particle density $`n_X(\text{r})`$ are

$$_d(E)=\frac{d𝒩(E,M_X/2)}{dE}\frac{1}{\tau _d}d^3\text{r}\frac{n_X(\text{r})}{2\pi |\text{r}_{}\text{r}|^2}$$
(1)
and

$$(E)=\frac{d𝒩(E,M_X)}{dE}\sigma _Avd^3\text{r}\frac{n_X(\text{r})^2}{2\pi |\text{r}_{}\text{r}|^2},$$
(2)
respectively. The shape of the spectrum is determined by the function $`d𝒩(E,E_{\text{jet}})/dE`$, which gives the number of sufficiently long-lived particles per energy interval emerging from a primary jet of energy $`E_{\text{jet}}`$. Birkel and Sarkar and Berezinsky and Kachelriess have used Monte Carlo simulations to calculate the spectrum from decay of superheavy relics . Fodor and Katz as well as Sarkar and Toldrà employed numerical integrations of the DGLAP evolution equations . The resulting spectrum in the interesting range between $`10^{19}`$eV and $`10^{20}`$eV is similar to the spectrum from the old modified leading log approximation (MLLA) , which was employed by Berezinsky et al. in their initial proposal of UHE cosmic rays from decay of superheavy dark matter decay .

We had also used the MLLA limiting spectrum in our proposal of wimpzilla annihilation  (Fig. 1).

The annihilation cross sections from the normalization in Fig. 1 were calculated for a smooth NFW  background halo. These cross sections would be in conflict with unitarity bounds on reaction cross sections for particles of mass $`M_X10^{12}`$GeV and therefore the flux from wimpzilla annihilation in the smooth halo component will (most likely) be negligible222There is also the possibility that wimpzillas form meta-stable bound states, similar to monopolonium. Then the contribution from the smooth halo component could be important and would resemble the anisotropy signature expected from monopolonium or superheavy dark matter decay.. However, a fraction around a few percent of the dark matter in the halo is expected to exist in dense isothermal dark matter sub-clumps333See e.g.  for evidence for sub-clumps from $`N`$-body simulations, and  for evidence from gravitational lensing., and due to large overdensities in the cores wimpzilla annihilation in these sub-clumps can comfortably generate the observed UHE flux with annihilation cross sections well below the unitarity bounds . This does not change the good fit of the annihilation model to the spectrum, which is a concequence of $`d𝒩(E,M_X)/dE`$, since the difference between events originating in the smooth and clumpy components affects only the overall factor $`d^3\text{r}n_X(\text{r})^2|\text{r}_{}\text{r}|^2`$ in (2).

Due to the negligibility of a possible contribution from the smooth halo component a clumpy source structure will be an unmistakable signature of collisional wimpzilla annihilation as a source of UHE cosmic rays. In particular, nearby sub-clumps will appear as hot spots in the sky (see Fig. 4 in Ref. ).

To calculate the expected overall anisotropy (after averaging over the clumps in the case of annihilation) one can integrate the local decay or annihilation rates entering in (1) and (2) along lines through the galactic halo, taking into account the galactic halo density and the expected sub-halo distribution, respectively.

The left panel in Fig. 2 shows the result for annihilation in a relatively broad subhalo distribution with the sub-clump density dropping off at a characteristic radius $`R_c^{cl}=20`$kpc (see  for a discussion of the expected sub-halo distribution).

Figure 2. The left panel shows the envelope of the UHE flux from annihilation in a galactic subhalo distribution as a function of the angle towards the galactic center, whereas the right panel shows the UHE flux from decay in a galactic NFW profile.

## 3 CONCLUSION

All top-down scenarios for generation of UHECRs from superheavy dark matter predict a large photonic component in cosmic ray primaries above $`E_{GZK}4\times 10^{19}`$eV, and a more or less dominant role of the galactic center:

Both annihilation in bound states (monopolonium or ”wimpzillonium” decay) and decay of superheavy relics imply that the dominant smooth dark matter halo distribution would be the primary source of UHECRs, whence the distribution of UHECR arrival directions should track the overall dark matter halo.

Collisional annihilation, on the other hand, would almost exclusively appear in the cores of local dark matter overdensities, and therefore would primarily track local overdensities. Here some domination of arrival directions pointing back towards the galactic center is expected because the sub-clump population is expected to drop off with distance from the galactic center. A nearby clump out of the line of sight to the galactic center could dominate the local signal and appear as a hot spot.

The southern site of the Auger observatory sees the galactic center, and therefore the prediction of a clumpy source distribution with some accumulation towards the galactic center should be testable after about three years of full operation of the southern site, which should yield of the order of 100 UHECR events with energies at and beyond $`10^{20}`$eV.


text:
# 1 Introduction

## 1 Introduction

Ising spin models have been used to describe a tremendous spectrum of physical and non-physical problems . Among others, they have proven useful to understand some aspect of granular media. In particular compaction dynamics under tapping .

Herein, we present a simple 1-$`d`$ Ising model to describe qualitative features of ripple formation. Previous experimental and theoretical works have underlined the primary role played by saltation in the emergence of sand ripples and dunes . Saltation means grain ejection and eolian transport over long distances. The hopping and rolling of grains have generally led to travelling ripple structures which in turn merge and grow.

In the following, an Ising spin description of the dynamics of ripples formation is presented. A 1-$`d`$ Hamiltonain is suggested. It is studied numerically with a Kawasaki dynamics. A mean-field treatment is also performed. A dynamical phase diagram is obtained. Results are discussed and compared with recent cellular automaton results  and with recent underwater experiments .

## 2 The model

A natural way to map a one-dimensional granular landscape into a spin model is to consider a 1-$`d`$ Ising model $`\{\sigma _i=\pm 1\}`$. Let the wind blowing from left to right. At each site $`i`$, $`\sigma _i=+1`$ if the local slope is positive (exposed to the wind). It is $`\sigma _i=1`$ if the local slope is negative (screened by the ripple crest and unaffected by the wind). Successive slopes being then associated to the succession of spin domains. Spin domains represent the ripple sides. A sketch of such a mapping is shown in Figure 1.

In Figure 1, a grain is first shifted over a distance $`\mathrm{}=7`$ by saltation. An avalanche then occurs to relaxe the surface. From the magnetic viewpoint, both extraction and relaxation of one grain produce two pairs of spin flips meaning that the whole dynamics goes by spin pairs. Moreover, it proves that the dynamics does not modify the total magnetization $`m=_i\sigma _i`$. To embody these physical effects charaterizing the saltation-avalanche mechanism, we consider the so-called Kawasaki dynamics . Having a 1-$`d`$ Ising system the equilibrium total magnetization is zero. In addition, it will stays zero since the Kawasaki dynamics preserves the total magnetization.

Avalanches which relaxe the granular surface can be viewed as small thermal fluctuations along the ripple sides. In our spin mapping such a scheme corresponds to a ferromagnetic nearest neigbor coupling $`J_{nn}`$ divided by $`k_BT`$ where $`T`$ is the temperature and $`k_B`$ the Boltzman’s constant. Moreover, saltation occurs only on the faces exposed to the wind and leads to a perturbation at the position $`i+\mathrm{}`$ where $`\mathrm{}`$ is the distance overwhich the sand is transported by saltation. As for avalanches, we mimic this effect by a ferromagnetic coupling $`J_s`$ which connects a site $`i`$ to a site $`i+\mathrm{}`$, divided by $`k_BT`$.

On this basis we consider the following dimensionless Hamiltonian

$$E=K_{nn}\underset{i=\mathrm{}}{\overset{+\mathrm{}}{}}\sigma _i\sigma _{i+1}K_s\underset{i=\mathrm{}}{\overset{+\mathrm{}}{}}\sigma _i\sigma _{i+\mathrm{}}\delta _{\sigma _i,+1},$$
(1)
where $`\delta _{x,y}`$ is the Kronecker function, $`K_{nn}\frac{J_{nn}}{k_BT}`$ and $`K_s\frac{J_s}{k_BT}`$. The Kronecker delta function in the second term assumes that only faces exposed to the wind (spins +1 for a wind blowing from left to right) play a role for saltation. These assumptions seem reasonable and not of drastic consequences at our present level of investigation.

## 3 Numerical simulations

Numerical simulations have been performed on a 1-$`d`$ system with periodic boundary conditions, i.e. a loop. The initial spin loop configuration is random with $`m=0`$. Thereafter, one positive spin and one negative spin are randomly selected at each time step. The energy variation $`\mathrm{\Delta }E`$ due to the flipping of both selected spins is calculated. Periodic boundary conditions are used in the energy calculation. The selected spins are then flipped with some probability $`p=\frac{\mathrm{exp}\mathrm{\Delta }E}{1+\mathrm{exp}\mathrm{\Delta }E}`$ as in classical Monte Carlo methods . The fact that two opposite spins are chosen at each time step preserves the $`m=0`$ value of the order parameter during the whole simulation, i.e. a Kawasaki dynamics.

Since the magnetization $`m`$ is zero and constant, a dynamical“order parameter” can be defined with the domain wall density $`\rho `$, i.e. the density of nearest neighboring antagonist pairs of spins. One should note that

$$\rho =\frac{\xi }{L}$$
(2)
where $`\xi `$ is the mean size of the spin domains. Loop lengths up to $`L=10000`$ spins have been used. Our numerical results has been checked against variation in the parameters ($`K_{nn}`$, $`K_s`$ and $`\mathrm{}`$). Finite-size effects as well as the spin pattern time evolution have also been analyzed. it is worth noting that $`\mathrm{}`$ and $`L`$ have been always chosen to be incommensurate in order to avoid any artefact due to the periodic boundary conditions.

Figure 2 presents typical space-time patterns in which domain walls are drawn for characterisitic lengths $`L=144`$ and $`\mathrm{}=21`$. Three different sets of parameter values $`K_{nn}`$ and $`K_s`$ are shown. Depending on the parameter values three different regimes are obtained. Either disorder (Figure 2a), or order with domains of roughly the size $`\mathrm{}/2`$ (Figure 2b), or also order with larger domains (Figure 2c). Domains having a size of the order $`\mathrm{}/2`$ are called thereby “ripples”. In parallel much larger domains are called “dunes” in reference to simulations in .

Figure 3 contains the domain wall density $`\rho `$ as a function of time for the three regimes of Figure 2. It is observed: (i) an erratic signal for the disordered regime presenting huge fluctuations of $`\rho `$, (ii) an apparently stable density $`\rho \mathrm{}/2`$ after a short transient regime and (iii) a low density regime ($`\rho 0`$) decorated with small fluctuations after some transient regime.

Figure 4 is the associated resulting dynamical phase diagram as function of $`(K_{nn},K_s)`$ obtained numerically for $`\mathrm{}>1`$. Three different dynamical regimes are thus found. A region of the disorder regime is represented in grey in the phase diagram. The boundary with the order regime is found to be stable whatever the values of $`\mathrm{}>1`$ and $`L`$. Curves $`\rho (K_{nn})`$ are shown in Figure 5 for different values of $`K_s`$ where one can clearly see a crossover around $`K_{nn}1`$. At this point, the domain size is equivalent to $`\mathrm{}/2`$. Moreover, this boundary is stable in time. This result could be astonishing since finite-size effects are expected. One should recall that $`m=0`$ corresponds to the equilibrium state of a one-dimensional Ising model and thus fluctuations and the domain wall evolution take place in the equilibrium state. Those are mainly driven by the Ising interaction $`K_{nn}`$.

For the order regime (white region in Figure 4), two distinct regions coexist: ripples (for $`K_s>>K_{nn}`$) and dunes (for $`K_{nn}>>K_s`$). The dashed line in Figure 4 represents the boundary between these regimes. However, this boundary moves slightly when the system size $`L`$ or the saltation length $`\mathrm{}`$ is changed. In fact, the ripple state is metastable. After very long times, ripples can suddenly merge into dunes (as seen by geologists). The larger $`K_s`$, the more stable ripples are.

Different initial correlations have been examined. Starting from a perfect antiferromagnetic situation, the system reaches the same regimes as described above. The dynamics is similar to what we have observed for the random inital situation. However, if one starts with a system containing only two large domains, i.e. a single giant dune, the system never produces ripples but remains unchanged in the order region of the dynamical phase diagram. In the disorder regime (in grey in Figure 4), the initial giant dune is however destroyed.

Finally, one should note the case $`\mathrm{}=1`$ is quite particular. In that case, spin domains are produced even when $`K_{nn}=0`$ for which a disorder regime is expected. As soon as $`\mathrm{}>1`$, the results summarized in Figure 4 are recovered.

## 4 Mean-field treatment

Eventhough the Hamiltonian (1) is rather simple, it does not allow for an easy exact solution. For instance, the classical method of the transfer matrix  leads to searching for $`2^{\mathrm{}+1}\times 2^{\mathrm{}+1}`$ matrices. In order to solve this problem, we have considered as a first approximation that the mean domain size $`\xi `$ is the relevant parameter. Fluctuations in $`\xi `$ are neglected.

Assuming an order with domains of size $`\xi `$, the first product $`\sigma _i\sigma _{i+1}`$ of Eq.(1) gives $`\frac{\xi 2}{\xi }`$ due to the presence of domain walls which are energetically defavorable. The second product $`\sigma _i\sigma _{i+\mathrm{}}`$ of Eq.(1) gives a periodically oscillating function of period $`2\xi `$. Thus, the energy per lattice site $`E/L`$ can be written as

$$\frac{E}{L}K_{nn}\left(\frac{\xi 2}{\xi }\right)\frac{K_s}{2}\mathrm{cos}\left(\frac{\pi \mathrm{}}{\xi }\right),$$
(3)
which is drawn in Figure 6 as a function of $`\xi `$ and for arbitrary values of the coupling energies $`K_{nn}`$ and $`K_s`$. Several minima are observed. The selection of various lengths $`\mathrm{}/2`$, $`\mathrm{}/3`$, $`\mathrm{}/4`$,… is thus expected from the above shape of $`E(\xi )`$. The more stable local minima is for $`\mathrm{}/2`$, i.e. it corresponds to ripples of size $`\mathrm{}`$.

For $`\xi >\mathrm{}/2`$, a large potential barrier separates the ripple state from the “giant dune” state. The asymmetry of this barrier implies the irreversibility of the “ripple $``$ dune” process. The ripple metastability is thus also obtained from a mean-field treatment in agreement with the above numerical results.

Nevertheless, the mean-field solution does not predict any disordered dynamical regime. Indeed, fluctuations of $`\xi `$ are neglected. Fluctuations destroy the order when $`K_s`$ is small.

## 5 Mapping onto a two-dimensional $`\mathrm{}\times \mathrm{}`$ static problem

Although ripple formation is a dynamical phenomenon, our study is generated via Monte Carlo (equilibrium-like) simulations. The dynamical aspect survives only because the relaxation towards equilibrium is slow and hampered by metastability.

Another remark concerns Eq.(1). Upon setting $`\delta _{\sigma _i,+1}=\frac{1}{2}(1\sigma _i)`$ in Eq.(1), multiplying out the $`\sigma `$’s, and using the fact that $`m=_i\sigma _i`$ one finds

$$E=K_{nn}\underset{i}{}\sigma _i\sigma _{i+1}K_s\underset{i}{}\sigma _i\sigma _{i+\mathrm{}}+\frac{1}{2}K_sm$$
(4)
in which the last term is an irrelevant constant that is set equal to zero later. Replace now each site index $`i`$ by two integer coordinates $`j`$ and $`k`$ uniquely defined by the conditions

$$ij\mathrm{}+k\text{ with }k\{0,1,\mathrm{},\mathrm{}1\}$$
(5)
and write $`\sigma _i=\sigma _{j,k}`$. In terms of these new variables the Hamiltonian Eq.(4) takes the form

$$E=K_{nn}\underset{j=\mathrm{}}{\overset{\mathrm{}}{}}\underset{k=0}{\overset{\mathrm{}1}{}}\sigma _{j,k}\sigma _{j,k+1}K_s\underset{j=\mathrm{}}{\overset{\mathrm{}}{}}\underset{k=0}{\overset{\mathrm{}1}{}}\sigma _{j,k}\sigma _{j+1,k}$$
(6)
provided we interpret $`(j,k)`$ as a site on an $`\mathrm{}\times \mathrm{}`$ cylindrical lattice with screw boundary conditions (CL/sbc) around the cylinder. Hence $`E`$ of Eq.(6) describes Onsager’s anisotropic nearest-neighbor Ising model on a CL/sbc, but with the additional constraint $`m=0`$.

This sheds an interesting light on the simulation results. The relevant facts are

(a) The constraint $`m=0`$ and the sbc are irrelevant for the thermodynamics of this model, which is therefore the same as that of the Onsager model without constraint and in magnetic field zero.

(b) The $`\mathrm{}\times \mathrm{}`$ lattice has no phase transition. But in the limit $`\mathrm{}\mathrm{}`$ this model does have the ferromagnetic Onsager phase transition. Of course for $`\mathrm{}`$ sufficiently large simulations will show an order/disorder transition-like behavior.

(c) The ground state $`𝒢_0`$, due to the condition $`m=0`$, necessarily has an interface. If the $`j`$ direction (the one parallel to the cylinder axis) is infinite, $`𝒢_0`$ is a cylinder which is positively magnetized for $`j>0`$ and negatively magnetized for $`j0`$. If the $`j`$ direction is finite (let’s say $`J<jJ`$ for some $`J`$), then $`𝒢_0`$ is still a low-lying state of energy $`E(𝒢_0)\mathrm{}K_s`$, but another low-lying state appears viz. the state $`𝒢_1`$ which has positive magnetization for $`\mathrm{}/2k\mathrm{}1`$ and negative magnetization for $`0k<\mathrm{}/2`$ ($`\mathrm{}`$ supposed even for convenience). Its energy is $`E(𝒢_1)2JK_{nn}`$. Both states are local energy minima in configuration space (because in both the interface length is at a local minimum), and, with the Kawasaki dynamics of this paper, the higher one of the two states is metastable with respect to the lower one. It is evident that $`𝒢_0`$ corresponds to the ”dune” and $`𝒢_1`$ to the ”rippled state”.

## 6 Connexion to recent underwater experiments

The physics underlying both ripple and dune formation is far from being understood making any analytic result valuable. Indeed, it has been believed that dune dynamics are of a different nature than ripple ones. Moreover the question of ripple growing to eventually produce a dune is still open.

Bagnold  distinguished ripples and ridges. Periodicity of the former ones was time independent while periodicity of the latter ones was growing with time. The present work suggests in a simple way that both patterns are two aspects of the same phenomenon since the ripple state seems to be metastable.

Recent experiments of underwater sand  have exhibited the formation of ripples. There, sudden merging of ripples do also occur in qualitative agreement with our results (see the upper part of the dynamical phase diagram of Figure 4). Moreover, these authors have found that a “transition” takes place between a flat landscape regime and ripple formation at a critical value of the Reynolds number Re. This feature is to be put in parallel to the transition that we have obtained (see Figure 4) between a disordered regime and a ripple dynamical regime. Indeed, the saltation length $`\mathrm{}`$ could be directly connected to the physical parameter Re.

## 7 Comparaison with a cellular automaton

In addition, a recent cellular automaton model  did produce a dynamical ripple state which eventually stabilizes into a giant dune. This model considers a granular landscape for which surface local fluctuations are dissipated through both saltation and avalanches . Saltation is restricted to grains exposed to the wind while avalanches occurs everywhere.

The associated dynamics is rather complex. First, ripples appear. Secondly, they coalesce very slowly into “giant dunes”. Main features of ripples and dunes formation are recovered within the present spin framework except for the disordered regime which is absent in .

Eventhough the cellular automaton model as well as actual spin-like model are far from real desert dunes and ripples, they seems to reproduce some analogous features.

## 8 Conclusion

In summary, we have proposed a 1-$`d`$ Ising model for ripple formation. The spin model has been studied numerically. In addition, a mean-field solution for the ripple size was also obtained. Moreover, the ripple state has been found to be metastable. A result which is in full agreement with very recent experiments under water as well as a previous landscape simulations.

Nature is obviously much more complex that a one-dimensional spin model. However our very simple description reproduces some basic features of dune formation as seen above. Moreover, the model could be extended to include, for instance, a distribution of saltation lengths or even, different grain types. This is left for future work. We do hope our results will stimulate additional experimental work in particular to check our dynamical phase diagram.

## 9 Acknowledgments

NV thanks the FNRS (Brussels, Belgium). He also thanks the LMDH for its hospitality during the progress of this work. We would like to thank Dr H. J. Hilhorst for a very stimulating comment.

Figure Captions

Figure 1 — (top) Sketch of the rule. A grain at position $`i`$ is unstable and is carried over by the wind to the position $`i+\mathrm{}`$ where an avalanche relaxes the surface. Spin signs indicate the slope of the landscape. (bottom) After the saltation-avalanche event, the magnetization remains unchanged.

Figure 2 — Typical space-time patterns in which domain walls are drawn: $`L=144`$ and $`\mathrm{}=21`$. Three different sets of parameter values are shown: (a) $`K_{nn}=0.5`$ and $`K_s=0.5`$, (b) $`K_{nn}=2.0`$ and $`K_s=2.5`$, and (c) $`K_{nn}=2.5`$ and $`K_s=0.5`$.

Figure 3 — The density of domain walls $`\rho `$ as a function of time. Three different sets of parameter values are shown, from top to bottom: $`K_{nn}=0.5`$ and $`K_s=0.5`$, $`K_{nn}=2.0`$ and $`K_s=2.5`$, and $`K_{nn}=2.5`$ and $`K_s=0.5`$.

Figure 4 — The dynamical phase diagram $`(K_{nn},K_s)`$ of the model. The disorder phase is represented in grey in the diagram. The order phase is represented in white. Ripples and dunes are also distinguished.

Figure 5 — The density of domain walls $`\rho `$ as a function of the Ising interaction $`K_{nn}`$. Three different values of $`K_s`$ are shown. The crossover separating disorder and order regimes is emphasized.

Figure 6 — Typical shape of the energy per lattice site $`E/L`$ as a function of the normalized domain size $`\xi /\mathrm{}`$.