text:
# Kondo effect in quantum dots coupled to ferromagnetic leads

## Abstract

We study the Kondo effect in a quantum dot which is coupled to ferromagnetic leads and analyse its properties as a function of the spin polarization of the leads. Based on a scaling approach we predict that for parallel alignment of the magnetizations in the leads the strong-coupling limit of the Kondo effect is reached at a finite value of the magnetic field. Using an equation-of-motion technique we study nonlinear transport through the dot. For parallel alignment the zero-bias anomaly may be split even in the absence of an external magnetic field. For antiparallel spin alignment and symmetric coupling, the peak is split only in the presence of a magnetic field, but shows a characteristic asymmetry in amplitude and position.

The Kondo effect hewson-book  in electron transport through a quantum dot (QD) with an odd number of electrons is experimentally well established goldhaber ; kondo-odd . Screening of the dot spin due to the exchange coupling with lead electrons yields, at low temperatures, a Kondo resonance. The main goal of the present work is to investigate how ferromagnetic leads influence the Kondo effect. In the extreme case of half-metallic leads, minority-spin electrons are completely absent, i.e., the screening of the dot spin is not possible, and no Kondo-correlated state can form. What happens, however, for the generic case of partially spin polarized leads? How does the spin-asymmetry affect the Kondo effect? Is there still a strong coupling limit, and how are transport properties modified?

Based on a poor man’s scaling analysis we first show that the strong-coupling limit can still be reached in this case if an external magnetic field is applied. This is familiar from the Kondo effect in QDs with an even number of electrons sasaki ; cobden ; kondo-even-theory-1 ; kondo-even-theory-2 , which occurs at finite magnetic fields, although the physical mechanism is different in the present case. In the second part of the paper we analyze within an equation-of-motion (EOM) approach the nonlinear transport through the QD. We find that for parallel alignment of the lead magnetizations the zero-bias anomaly is split. This splitting can be removed by appropriately tuning the strength of an external magnetic field $`B`$. In the antiparallel configuration of the lead magnetizations no splitting occurs at zero field.

The Anderson Hamiltonian for a QD with a single level at energy $`ϵ_0`$ coupled to ferromagnetic leads is

$`H`$ $`=`$ $`{\displaystyle \underset{rk\sigma }{}}\epsilon _{rk\sigma }c_{rk\sigma }^{}c_{rk\sigma }+ϵ_0{\displaystyle \underset{\sigma }{}}d_\sigma ^{}d_\sigma +Ud_{}^{}d_{}d_{}^{}d_{}`$ (1)
$`+{\displaystyle \underset{rk\sigma }{}}(V_{rk}d_\sigma ^{}c_{rk\sigma }+V_{rk}^{}c_{k\sigma }^{}d_\sigma )+g\mu _BBS_z,`$
where $`c_{rk\sigma }`$ and $`d_\sigma `$ are the Fermi operators for electrons with wavevector $`k`$ and spin $`\sigma `$ in the leads, $`r=L,R`$, and in the QD, $`V_{rk}`$ is the tunneling amplitude, $`S_z=(d_{}^{}d_{}d_{}^{}d_{})/2`$, and the last term is the Zeeman energy of the dot. (Stray fields from the leads are neglected.) We assume identical leads and symmetric coupling, $`V_{Lk}=V_{Rk}`$. The ferromagnetism of the leads is accounted for by different densities of states (DOS) $`\nu _r(\omega )`$ and $`\nu _r(\omega )`$ for up and down-spin electrons.

In the following we study the two cases of parallel (P) and antiparallel (AP) alignment of the leads’ magnetic moments. For the AP configuration and zero magnetic field and bias voltage, the model is equivalent (by canonical transformation glazman ) to a QD coupled to a single lead with DOS $`\nu _L+\nu _R=\nu _L+\nu _R`$. In this case, the usual Kondo resonance forms, which is the same as for nonmagnetic electrodes hewson-book .

This changes for the P configuration. In this case, there is an overall asymmetry for up and down spins, say $`\nu _L+\nu _R>\nu _L+\nu _R`$. To understand how this asymmetry affects the Kondo physics we apply the poor man’s scaling technique anderson , performed in two stages haldane . In the first stage, charge fluctuations dominate and lead to a renormalization of the QD’s levels. Since the renormalization for the spin-down level is stronger than for spin-up, a level splitting between the two spin orientations is generated. This is one of the key mechanism for all the effects discussed below. To reach the strong-coupling limit it is, therefore, essential to apply an external magnetic field to compensate for the generated spin splitting. In the second stage, the resulting model is mapped onto a Kondo Hamiltonian and the degrees of freedom involving spin fluctuations are integrated out. For simplicity we assume for the scaling analysis flat DOSs $`\nu _{r\sigma }`$ and neglect the $`k`$-dependence of the tunnel amplitudes $`V_{rk}=V`$.

First we reduce the cutoff $`D`$ from $`D_0`$, which is the smaller value of the bandwidth and the onsite repulsion $`U`$ haldane . Charge fluctuations lead to the scaling equations

$`{\displaystyle \frac{dϵ_\sigma }{d\mathrm{ln}(D_0/D)}}={\displaystyle \frac{\mathrm{\Gamma }_{\overline{\sigma }}}{2\pi }}`$ (2)
where we defined $`\mathrm{\Gamma }_\sigma =2\pi |V|^2_r\nu _{r\sigma }`$, and $`\overline{\sigma }`$ is opposite to $`\sigma `$. This yields the solution $`\mathrm{\Delta }\stackrel{~}{ϵ}=\stackrel{~}{ϵ}_{}\stackrel{~}{ϵ}_{}=(1/\pi )P\mathrm{\Gamma }\mathrm{ln}(D_0/D)+\mathrm{\Delta }ϵ_0`$, where $`P=(\nu _r\nu _r)/(\nu _r+\nu _{r)}`$ measures the spin polarization in the leads, $`\mathrm{\Gamma }=(\mathrm{\Gamma }_{}+\mathrm{\Gamma }_{})/2`$ and $`\mathrm{\Delta }ϵ_0=g\mu _BB`$ is the Zeeman splitting. The empty-dot state $`|0`$ hybridizes with states where the dot is singly occupied $`|1\sigma `$ with either spin up or down, while the singly-occupied state $`|1\sigma `$ only hybridizes with the empty-dot state $`|0`$ (for $`U|ϵ|`$ \- asymmetric Anderson model). Because of the spin-dependent DOS in the leads the hybridization is spin-dependent, which is the physical origin of the generated $`\mathrm{\Delta }\stackrel{~}{ϵ}`$.

To describe Kondo physics (for $`\stackrel{~}{ϵ}<0`$) we terminate haldane  the scaling of Eq. (2) at $`\stackrel{~}{D}\stackrel{~}{ϵ}`$, and perform a Schrieffer-Wolff transformation. Using the renormalized parameters $`\stackrel{~}{D}`$ and $`\stackrel{~}{ϵ}`$ we get the effective Kondo Hamiltonian

$`H_{\mathrm{Kondo}}={\displaystyle \underset{rr^{}kk^{}}{}}\{J_+S^+c_{rk}^{}c_{r^{}k^{}}+J_{}S^{}c_{rk}^{}c_{r^{}k^{}}+`$
$`S_z(J_zc_{rk}^{}c_{r^{}k^{}}J_zc_{rk}^{}c_{r^{}k^{}})\},`$ (3)
plus the term $`2S_z\stackrel{~}{D}(\nu _{}J_z\nu _{}J_z)`$ and a potential scattering term. The initial values for the coupling constants are $`J_+=J_{}=J_z=J_z=|V|^2/|\stackrel{~}{ϵ}|J_0`$. To reach the strong-coupling limit we tune the external magnetic field $`B`$ such that the total effective Zeeman splitting vanishes, $`\mathrm{\Delta }\stackrel{~}{ϵ}=0`$ (the field $`B`$ will also slightly modify the DOS in the leads kondo-even-theory-1 ). During the second stage of scaling spin fluctuations will renormalize the three coupling constants $`J_+=J_{}J_\pm `$, $`J_z`$, and $`J_z`$ differently. The scaling equations are

$`{\displaystyle \frac{d(\nu _\pm J_\pm )}{d\mathrm{ln}(\stackrel{~}{D}/D)}}`$ $`=`$ $`\nu _\pm J_\pm (\nu _{}J_z+\nu _{}J_z)`$ (4)
$`{\displaystyle \frac{d(\nu _\sigma J_{z\sigma })}{d\mathrm{ln}(\stackrel{~}{D}/D)}}`$ $`=`$ $`2(\nu _\pm J_\pm )^2`$ (5)
with $`\nu _\pm =\sqrt{\nu _{}\nu _{}}`$, $`\nu _\sigma _r\nu _{r\sigma }`$ com\_models . To solve these equations we observe that $`(\nu _\pm J_\pm )^2(\nu _{}J_z)(\nu _{}J_z)=0`$ and $`\nu _{}J_z\nu _{}J_z=J_0P(\nu _{}+\nu _{})`$ is constant as well. I.e., there is only one independent scaling equation. All coupling constants reach the stable strong-coupling fixed point $`J_\pm =J_z=J_z=\mathrm{}`$ at the Kondo energy scale, $`Dk_BT_K`$. For the P configuration the Kondo temperature in leading order,

$$T_\mathrm{K}(P)\stackrel{~}{D}\mathrm{exp}\left\{\frac{1}{(\nu _{}+\nu _{})J_0}\frac{\mathrm{arctanh}(P)}{P}\right\},$$
(6)
depends on the polarization $`P`$ in the leads. It is a maximum for nonmagnetic leads, $`P=0`$, and vanishes for $`P1`$.

Finally, we point out an interesting consequence of the spin polarization in the leads. With nonmagnetic leads, the Kondo Hamiltonian couples the spin of the QD to the spin of the leads only, but not to its charge. To analyze the analogous situation in our case, we introduce the (pseudo) spin $`\stackrel{}{\sigma }=(1/2)_{kk^{}rr^{}\sigma \sigma ^{}}c_{kr\sigma }^{}\sigma _{\sigma \sigma ^{}}c_{k^{}r^{}\sigma ^{}}/(2\stackrel{~}{D}\sqrt{\nu _\sigma \nu _\sigma ^{}})`$, where the spin-dependent normalization factor is crucial to ensure the proper spin commutation relations, and the (pseudo) charge $`en=e_{kk^{}rr^{}\sigma }c_{kr\sigma }^{}c_{k^{}r^{}\sigma }/(2\stackrel{~}{D}\nu _\sigma )`$. The last term in Eq. (3) can, then, be written as $`2\stackrel{~}{D}(\nu _{}J_z+\nu _{}J_z)\sigma _zS_z`$ plus $`\stackrel{~}{D}(\nu _{}J_z\nu _{}J_z)nS_z`$. The first term is analogous to the Kondo model with nonmagnetic leads, while the second term couples spin to charge. The latter does not scale up and the associated additional renormalization of the Zeeman splitting, $`(1/\pi )P\mathrm{\Gamma }`$, is negligible as compared to $`\mathrm{\Delta }\stackrel{~}{ϵ}`$ in the limit $`D_0|ϵ|`$.

The unitary limit for the P configuration can be achieved by tuning the magnetic field appropriately, as discussed above. In this case, the maximum conductance through the QD is $`G_{\mathrm{max},\sigma }^P=e^2/h`$ per spin, i.e., the same as for nonmagnetic leads. This yields that the amplitude of the Kondo resonance for up and down spin at the Fermi level are different, since $`G_{\mathrm{max},\sigma }^P\mathrm{\Gamma }_\sigma (E_\mathrm{F})\rho _\sigma (E_\mathrm{F})`$, and therefore, $`\rho _{}(E_\mathrm{F})/\rho _{}(E_\mathrm{F})=(1P)/(1+P)`$. For the AP configuration, the maximal conductance is reduced, $`G_{\mathrm{max},\sigma }^{AP}=(1P^2)e^2/h`$, and vanishes for $`P1`$.

In the remainder of this article we analyze the DOS of the QD and address nonequilibrium transport. For a qualitative discussion we should employ the simplest technique which accounts for both the formation of Kondo resonances and the influence of the spin-dependent renormalization of the dot level on spin fluctuations. The equation-of-motion (EOM) technique with the usual decoupling procedure meir ; EOM  for higher order Green functions (GF) satisfies the first requirement but not the second. We, therefore, extend this scheme by calculating the level splitting $`\mathrm{\Delta }\stackrel{~}{ϵ}`$ self-consistently. For a more quantitative analysis one could include higher-order (than usually) Green’s functions in the EOM approach or higher-order diagrams in the resonant-tunneling-approximation konig , or use more advanced schemes such as real-time konig2  or numerical RG costi1  methods. These techniques are, however, much more complex com\_slaveboson  and are not pursued here.

Within the Keldysh formalism, the transport current $`I=_\sigma I_\sigma `$ through a QD for $`\mathrm{\Gamma }_{\mathrm{R}\sigma }(\omega )=\lambda _\sigma \mathrm{\Gamma }_{\mathrm{L}\sigma }(\omega )`$ is

$`I_\sigma ={\displaystyle \frac{e}{\mathrm{}}}{\displaystyle 𝑑\omega \frac{\mathrm{\Gamma }_{\mathrm{L}\sigma }(\omega )\mathrm{\Gamma }_{\mathrm{R}\sigma }(\omega )}{\mathrm{\Gamma }_{\mathrm{L}\sigma }(\omega )+\mathrm{\Gamma }_{\mathrm{R}\sigma }(\omega )}[f_\mathrm{L}(\omega )f_\mathrm{R}(\omega )]\rho _\sigma (\omega )}`$ (7)
where $`\rho _\sigma (\omega )=(1/\pi )\mathrm{Im}G_\sigma ^{\mathrm{ret}}(\omega )`$. For strong interaction ($`U\mathrm{}`$) the retarded Green’s function is

$`G_\sigma ^{\mathrm{ret}}(\omega )={\displaystyle \frac{1n_{\overline{\sigma }}}{\omega ϵ_\sigma \mathrm{\Sigma }_{0\sigma }(\omega )\mathrm{\Sigma }_{1\sigma }(\omega )+i0^+}}`$ (8)
where $`\mathrm{\Sigma }_{0\sigma }(\omega )=_{k\mathrm{L},\mathrm{R}}|V_k|^2/(\omega \epsilon _{k\sigma })`$ is the self-energy for a noninteracting QD, while

$`\mathrm{\Sigma }_{1\sigma }(\omega ,\mathrm{\Delta }\stackrel{~}{ϵ})={\displaystyle \underset{k\mathrm{L},\mathrm{R}}{}}{\displaystyle \frac{|V_k|^2f_{\mathrm{L}/\mathrm{R}}(\epsilon _{k\overline{\sigma }})}{\omega \sigma \mathrm{\Delta }\stackrel{~}{ϵ}\epsilon _{k\overline{\sigma }}+i\mathrm{}/2\tau _{\overline{\sigma }}}}`$ (9)
appears for interacting QDs only. The average occupation of the QD with spin $`\sigma `$ is obtained from $`n_\sigma =(i/2\pi )𝑑\omega G_\sigma ^<(\omega )`$. Extending the standard derivation meir , we replaced on the r.h.s. of Eq. (9) $`\mathrm{\Delta }ϵ\mathrm{\Delta }\stackrel{~}{ϵ}`$, where $`\stackrel{~}{ϵ}_\sigma `$ is found self-consistently from the relation

$`\stackrel{~}{ϵ}_\sigma =ϵ_\sigma +\mathrm{Re}[\mathrm{\Sigma }_{0\sigma }(\stackrel{~}{ϵ}_\sigma )+\mathrm{\Sigma }_{1\sigma }(\stackrel{~}{ϵ}_\sigma ,\mathrm{\Delta }\stackrel{~}{ϵ})],`$ (10)
which describes the renormalized dot-level energy, where the real part of the denominator of Eq. (8) vanishes hewson-book . We emphasize that without this self-consistency relation the Kondo resonances will, in general, appear at different positions, which disagree with the conclusions from the scaling analysis sergueev . The procedure simulates higher-order contributions and the influence of the renormalization of the dot level on spin fluctuations. Following Ref. meir  we introduce, in heuristic way, a lifetime $`\tau _\sigma (\mu _\mathrm{L},\mu _\mathrm{R},\stackrel{~}{ϵ}_{},\stackrel{~}{ϵ}_{})`$ which describes decoherence due to a finite bias voltage $`V`$ or level splitting $`\mathrm{\Delta }\stackrel{~}{ϵ}`$. It is obtained in second-order perturbation theory and depends on the electrochemical potentials in the leads, $`\mu _{\mathrm{L}(\mathrm{R})}`$, and the level positions. Again, we replace the bare levels by the renormalized ones. In the numerical results presented below we use Lorentzian bands of width $`D=100\mathrm{\Gamma }`$.

For nonmagnetic leads, $`P=0`$ and zero magnetic field, $`B=0`$, the proposed approximation is identical to the standard EOM scheme meir . For finite magnetic field, $`B0`$, the self-consistency condition yields a splitting of the Kondo resonances which is slightly smaller than $`2g\mu _BB`$, in agreement with both experimental goldhaber  and theoretical findings moore ; costi2 . For $`B=0`$ and $`P>0`$ in the parallel configuration, we obtain a value of the splitting $`\mathrm{\Delta }\stackrel{~}{ϵ}`$ comparable to the result from scaling, Eq. (2).

In Fig. 1 we plot the DOS of the QD for spins in AP and P configurations with spin polarization $`P=0.2`$ in the leads. In the AP configuration there is one Kondo resonance \[Fig. 1(a)\] and the DOS is the same as for the case of nonmagnetic leads. For the P configuration, however, the Kondo resonance splits \[Fig. 1(c)\], which can be compensated by an external magnetic field $`B`$ \[Fig. 1(d)\]. In the latter case, the amplitude of the Kondo resonance for spin down significantly exceeds that for spin up (as discussed above). A finite bias voltage, $`\mu _\mathrm{R}\mu _\mathrm{L}=eV>0`$, again leads to a splitting for both the AP and the P configurations, \[Fig. 1(b) and (e)\]. In the AP configuration, the amplitude of the upper and the lower Kondo peak appear asymmetric \[Fig. 1(b)\].

In Fig. 2 we show the differential conductance as a function of the transport voltage. For nonmagnetic leads, there is a pronounced zero-bias maximum \[Fig. 2(a)\], which splits in the presence of a magnetic field \[Fig. 2(b)\]. For magnetic leads and parallel alignment, we find a splitting of the peak in the absence of a magnetic field \[Fig. 2(c)\], which can be tuned away by an appropriate magnetic field \[Fig. 2(d)\]. In the AP configuration, the opposite happens, no splitting at $`B=0`$ \[Fig. 2(e)\] but finite splitting at $`B>0`$ \[Fig. 2(f)\] with an additional asymmetry in the peak amplitudes as function of the bias voltage. This asymmetry is related to the asymmetry in the amplitude of DOS \[Fig. 1(b)\]. All these findings are in good agreement with our scaling analysis. In Fig. 2(g) we show the tunnel magnetoresistance (TMR) which can be much larger than for conventional TMR systems. Finally, we find that the positions of the peaks in the AP configuration in the presence of a magnetic field are slightly shifted as a function of the polarization $`P`$ \[Fig. 2(h)\]. This can be explained in the similar way as in Ref. kondo-even-theory-2  by an additional level splitting $`\delta \mathrm{\Delta }\stackrel{~}{ϵ}=(1/4\pi )P\mathrm{\Gamma }/\stackrel{~}{ϵ}eV`$ at finite bias voltages due to spin accumulation in the QD.

We finally comment on the observability of our proposal. Ferromagnetic systems usually have strong spin-orbit coupling which could lead to spin-flip relaxation and, thus, could destroy the Kondo effect. On the other hand, recent observation of Aharonov-Bohm oscillations in a ferromagnetic ring kasai  proved coherent transport, i.e., the formation of the Kondo cloud in ferromagnetic leads might be possible as well. How can one attach ferromagnetic leads to a QD? A conceivable realization might be to put carbon nanotubes in contact to ferromagnetic leads ago . The Kondo effect has been observed already for nonmagnetic electrodes cobden ; bachtold . Alternatively one might use magnetic tunnel junctions with magnetic impurities in the barrier, or spin-polarized STM ralph ; wiesendanger .

In conclusion, we presented a qualitative study of the Kondo effect in QDs coupled to ferromagnetic leads. In particular we found a splitting of the Kondo resonance for parallel alignment of the leads magnetizations, even in the absence of a magnetic field. Our results are based on a poor man’s scaling approach and an EOM technique. Further investigations on a more quantitative level using more advanced techniques would be desirable.

We thank L. Borda, J. von Delft, L. Glazman, B. Jones, Yu.V. Nazarov, A. Rosch, H. Schoeller, A. Tagliacozzo, M. Vojta and A.Zawadowski for discussion. Support by the Research Project KBN 5 P03B 091 20 and the DFG through the CFN and the Emmy-Noether program is acknowledged.


format: text

text:
# Landau Expansion for the Critical Point of a Polydisperse System

## I Introduction

Many systems, such as colloidal suspensions, and polymer blends and solutions, contain particles with an effectively continuous distributions of an attribute ($`\sigma `$, say) such as size, charge, or chemical character. Such systems are termed polydisperse . Their phase equilibria can become complicated, but may exhibit relatively simple features including coexistence of two liquid phases of different composition. In some cases the tendency to phase separate is inherited from the underlying monodisperse physics (e.g. attractive colloids); in others it stems only from polydispersity (e.g. fractionation among chains with a spectrum of chemical compositions).

In a polydisperse system undergoing fluid-fluid demixing, the phase diagram looks rather different from that of a simple binary mixture even if the underlying physics is similar. It is common experimentally to draw the phase diagram on a composition/temperature plane; one effect of polydispersity is then to split the binodal curve into cloud and shadow curves that cross at the critical point. (The latter is not at the maximum of either curve.) This is because, when a system is cooled to the point where it becomes globally unstable to phase separation (the cloud point), the incipient new phase (its shadow) differs from the old one not only in the number of particles per unit volume, as usual, but in its $`\sigma `$-distribution. Because of this difference, the incipient new phase does not itself represent a point on the cloud curve. (However, it does lie on the cloud curve of a different system whose overall size distribution coincides with that of the incipient, shadow phase).

An interesting question is whether all this behavior can be captured within a Landau expansion. There are at least two reasons to consider this. Firstly, conceptual (as well as practical) progress has been made recently in understanding phase equilibria in polydisperse systems within a certain class of mean-field models: those in which the excess free energy depends on the $`\sigma `$ distribution only through a finite set of moments . A promising method involves projecting the free energy onto a relatively small number of densities; when these are chosen correctly, the cloud, shadow, and spinodal curves, and with them any critical points, are found without approximation beyond those used to set up the mean-field free energy in the first place.

Due to their mean-field character and restricted set of relevant densities, such projection techniques should be reconcilable with a Landau expansion near the critical point – yet the unusual topology of the phase diagram appears, at first sight, to preclude this. In particular, the usual adoption of a quartic expansion in one density (with coefficients dependent only on temperature) is bound to exhibit a critical point at the apex of a binodal curve. We show below that a more careful Landau expansion does not suffer from this restriction, and instead reproduces the correct, more complicated behavior outlined above.

A second reason to study the Landau expansion is in the hope of going beyond mean field theory. Experimental values for the critical exponents of a polydisperse system  deviate from those of the monodisperse case and are consistent with Fisher-renormalised exponents . However, to our knowledge, there is no explicit renormalization group calculation for a polydisperse system. Also, quite apart from calculating critical indices, it is important for many purposes (see e.g. ) to properly classify and identify the scaling fields near the critical point, and this task cannot really be begun without a consistent description at the Landau level. Hence the Landau expansion presented below represents a first step towards a fuller understanding of polydisperse critical behavior beyond mean field level.

## II Monodisperse system

In this section we review briefly the classical Landau expansion of the free energy around the critical point for a simple monodisperse system. This will facilitate a comparison with the polydisperse case and establish the notation.

In a canonical ensemble, the free energy $``$ is defined as a function of the temperature $`T`$, the number of particles $`N`$ and the volume $`V`$. In the thermodynamic limit, however, we can define a reduced free energy (density) as

$`F(\rho ,T)={\displaystyle \frac{(N,T,V)}{Vk_BT}},`$ (1)
which is a function only of the density of particles $`\rho =N/V`$ and the temperature. At fixed temperature, any terms that are constant or linear in density are irrelevant for obtaining the phase equilibrium of the system (which is our main aim) since they can be assimilated into shifts of the entropy and of the chemical potential.

### A Critical Point

Within the mean field approximation, which is where the Landau expansion holds, the critical point must obey one condition common to all points on the spinodal curve,

$`{\displaystyle \frac{^2F}{\rho ^2}}|_c=0,`$ (2)
and a second condition particular to the critical point,

$`{\displaystyle \frac{^3F}{\rho ^3}}|_c=0,`$ (3)
where $`|_c`$ means evaluated at the critical point $`(\rho ^c,T_c)`$.

### B Landau Expansion

Once the critical point is determined, the expansion follows as usual:

$`F(\rho ,T){\displaystyle \frac{1}{2}}{\displaystyle \frac{^2F}{\rho ^2}}|_{\rho ^c,T}\delta \rho ^2+`$ (4)
$`{\displaystyle \frac{1}{3!}}{\displaystyle \frac{^3F}{\rho ^3}}|_{\rho ^c,T}\delta \rho ^3+{\displaystyle \frac{1}{4!}}{\displaystyle \frac{^4F}{\rho ^4}}|_{\rho ^c,T}\delta \rho ^4+\mathrm{},`$ (5)
where $`\delta \rho \rho \rho ^c`$. We forget about the constant and linear terms since, as described above, they can be viewed as (temperature dependent) shifts in the entropy and chemical potential with no effect on phase behavior.

The coefficients of the expansion depend on the temperature and can be expanded themselves about the critical temperature. Using Eqs. (2) and (3), and keeping the lowest order term for each coefficient, we obtain:

$`F(\rho ,T){\displaystyle \frac{1}{2}}{\displaystyle \frac{^3F}{\rho ^2T}}|_c\delta T\delta \rho ^2+`$ (6)
$`{\displaystyle \frac{1}{3!}}{\displaystyle \frac{^4F}{\rho ^3T}}|_c\delta T\delta \rho ^3+{\displaystyle \frac{1}{4!}}{\displaystyle \frac{^4F}{\rho ^4}}|_c\delta \rho ^4+\mathrm{},`$ (7)
where $`\delta TTT_c`$. Since $`\delta T\delta \rho ^2`$ (see Eq. (9) below) the term in $`\delta T\delta \rho ^3`$ is of higher order and can be neglected. Therefore, the lowest order expansion of the free energy has the simple form:

$`F(\rho ,T)A\delta T\delta \rho ^2+C\delta \rho ^4,`$ (8)
for appropriate $`A`$ and $`C`$. The coexisting phases obey the standard result

$`\delta T{\displaystyle \frac{2C}{A}}\delta \rho ^2,`$ (9)
giving the well-known parabolic shape of the phase diagram around the critical point for monodisperse systems (always within the mean field approximation).

## III Polydisperse System

In the polydisperse case, the particles are not all identical but differ in certain scalar attributes (particle radius, charge, chain length…) which can be tagged with indices $`\sigma _1,\sigma _2,\mathrm{}`$. For simplicity, we focus on systems with just one polydisperse index $`\sigma `$, and write the corresponding density of particles as $`\rho (\sigma )`$. In this case, the free energy can be written as two parts, an ideal term and an excess term, as follows :

$`F[\rho ]={\displaystyle 𝑑\sigma \rho (\sigma )\left\{\mathrm{log}\left(\rho (\sigma )\right)1\right\}}+F_{ex}[\rho ].`$ (10)
While the ideal part only depends on $`\rho (\sigma )`$ as shown, the excess part $`F_{ex}`$ is, in general, an arbitrary functional of $`\rho (\sigma )`$ and a function of $`T`$. Here, in order to explore how polydispersity affects the Landau expansion around the critical point in a tractable manner, we restrict ourselves to cases where the functional dependence of the excess part of the free energy on the density distribution $`\rho (\sigma )`$ is through a single ‘generalized moment’  of the distribution:

$`F[\rho ]={\displaystyle 𝑑\sigma \rho (\sigma )\left\{\mathrm{log}\left(\rho (\sigma )\right)1\right\}}+f(\rho _1,T),`$ (11)
where

$`\rho _1={\displaystyle 𝑑\sigma \rho (\sigma )w_1(\sigma )}`$ (12)
with $`w_1(\sigma )`$ an arbitrary weight function. Although this assumption seems restrictive, it is true for some important models including the Flory-Huggins theory of length-polydisperse polymer solutions and blends . The models for which the excess part of the free energy depends on the density $`\rho (\sigma )`$ through a finite set of moments of this distribution, for which Eq. (11) is a particular example, are usually termed truncatable models.

### A The Parent Distribution Function

Polydispersity typically appears in the process of chemical manufacture of the particles of the system, so that the relative proportion of particles with different $`\sigma `$ is fixed thereafter. The actual density of each species then depends on how much solvent is present. We define a parent distribution $`P(\sigma )`$ as

$`P(\sigma )=P_0n(\sigma ),`$ (13)
where $`n(\sigma )`$ is a normalized distribution (fixed by the initial preparation of the system) and $`P_0`$ is the overall density of particles, which varies as one moves along a ‘dilution line’ through the phase diagram. All other moments of the parent are fixed by $`P_0`$; for example,

$`P_1={\displaystyle 𝑑\sigma P(\sigma )w_1(\sigma )}=P_0n_1`$ (14)
where $`n_1`$ is fixed at synthesis via the parent:

$`n_1={\displaystyle 𝑑\sigma n(\sigma )w_1(\sigma )}.`$ (15)
Although the parent distribution is globally fixed, the system can separate into coexisting phases with a distribution $`\rho (\sigma )`$ that is locally different from $`P(\sigma )`$. If there are two such phases $`\alpha `$ and $`\beta `$, with respective distributions $`\rho ^\alpha (\sigma )`$ and $`\rho ^\beta (\sigma )`$, the parent distribution is recovered as the volumetric average:

$`x^\alpha \rho ^\alpha (\sigma )+x^\beta \rho ^\beta (\sigma )=P(\sigma ),`$ (16)
where $`Vx^\alpha `$ and $`Vx^\beta `$ are the volumes of phases $`\alpha `$ and $`\beta `$, respectively. Note however that at a phase boundary which marks the onset of phase separation, the volume fraction of the incipient (shadow) phase is negligible and, therefore, the $`\sigma `$-distribution of the original (cloud) phase must coincide with the parent distribution.

### B Projection Method

Although for a given $`n(\sigma )`$ we require only one density (such as $`\rho _1`$), alongside temperature, to specify a point on the phase diagram, the corresponding free energy must depend on the full density distribution $`\rho (\sigma )=P_0n(\sigma )`$. Accordingly it is not correct to substitute the parent distribution in a functional like Eq. (11) to obtain a one-variable free energy $`F(P_0)`$ (or equivalently, $`F(\rho _1)`$). In fact such a substitution would describe a case where the normalized distribution $`n(\sigma )`$ were constrained to be equal in all coexisting phases – which is not the physical situation. Instead, in the case of Eq. (11), the equality of chemical potentials dictates that the $`\sigma `$-distribution of any cloud phase (which coexists, by definition, with an unperturbed parent) must obey

$`\rho (\sigma )=P(\sigma )e^{\lambda w_1(\sigma )},`$ (17)
for an appropriate value of $`\lambda `$. This distribution is not proportional to the parent unless $`\lambda =0`$ (when the phases are identical).

However, under certain conditions there is a procedure to transform the free energy functional into an effective free energy function. The conditions are first that (as in Eq. (11)) the excess free energy $`f`$ only depends on the distribution $`\rho (\sigma )`$ through a finite set of moments (defined as in Eq. (12)); and second that we are only interested in phase boundaries (the cloud and the shadow curves) and the spinodals. The procedure, called the projection method , involves several steps briefly described here, and leads to a free energy function dependent only on the moments required to specify the excess free energy.

First, considers the free energy

$`F[\rho ]={\displaystyle 𝑑\sigma \rho (\sigma )\left\{\mathrm{log}\left(\frac{\rho (\sigma )}{P(\sigma )}\right)1\right\}}+f(\rho _1,T),`$ (18)
which, for our purposes, is equivalent to Eq. (11), since the inclusion of the parent as a denominator in the logarithm represents a linear term (which can be considered as a mere shift of the chemical potentials). But note that, for that equivalence to hold, the parent $`P(\sigma )`$, including $`P_0`$, must be kept fixed as the density distribution $`\rho (\sigma )`$ varies: this observation is crucial to the Landau expansion developed below.

The projection method follows by minimising the first (ideal) term of this free energy subject to the constraint of fixed values for the moment(s) upon which the excess free energy depends. (The only such moment is $`\rho _1`$ in our case.) This yields

$`\rho (\sigma )=P(\sigma )e^{\lambda _1w_1(\sigma )},`$ (19)
where $`\lambda _1`$ is the Lagrange multiplier used to constrain the minimization, determined by $`\rho _1`$ through

$`\rho _1={\displaystyle 𝑑\sigma P(\sigma )w_1(\sigma )e^{\lambda _1w_1(\sigma )}}.`$ (20)
Since both the parent phase (the cloud) and the coexisting phase (the shadow, Eq. (17)) obey Eq. (19) for different values of $`\lambda _1`$, the method contains no approximation when it comes to finding the phase boundaries. However, if the $`\sigma `$-distribution of any of the phases adopts a form different from (19), the method produces only approximate results. This occurs whenever the system separates into finite volumes of more than one phase, that is, everywhere within the cloud curve. (In Eq. (16), the projection method is inexact when $`x_\alpha `$ and $`x_\beta `$ are both nonzero.)

The free energy obtained by this method is

$`F(\rho _1,T)=\lambda _1\rho _1\rho _0+f(\rho _1,T),`$ (21)
where

$`\rho _0={\displaystyle 𝑑\sigma P(\sigma )e^{\lambda _1w_1(\sigma )}}.`$ (22)
This free energy depends on only one density variable, since both $`\lambda _1`$ and $`\rho _0`$ are given by $`\rho _1`$. The coexisting phases are obtained from this free energy function via the common tangent method, as in the monodisperse case.

### C Dilution Line

As mentioned above, the projected free energy requires knowledge of the parent distribution function $`P(\sigma )`$, whose normalization $`P_0`$ varies as one moves along a dilution line on the phase diagram. For this reason, the projected free energy $`F(\rho _1,T)`$ depends parametrically on $`P_0`$, through $`\lambda _1(\rho _1,P_0)`$ and $`\rho _0(\lambda _1,P_0)`$; see Eqs. (20) and (22), respectively. This dependence is parametric because, in the process of calculating the coexistence between phases, $`P_0`$ must remain fixed, so as to maintain the equivalence of Eq. (11) and (18).

Suppose that, for fixed $`P_0`$ and $`T`$, we obtain by this method a coexistence between phases $`\alpha `$ and $`\beta `$, with moment densities $`\rho _1^\alpha `$ and $`\rho _1^\beta `$, respectively. This determines the value of $`\lambda _1^\alpha `$, $`\lambda _1^\beta `$, $`\rho _0^\alpha `$, and $`\rho _0^\beta `$ through Eqs. (20) and (22) which, in turn, determines the density distribution of the two phases $`\rho ^\alpha (\sigma )`$ and $`\rho ^\beta (\sigma )`$ via Eq. (19). In general, neither of these distributions is the parent distribution.

This is consistent with the fact that the parent distribution does not coexist with a second phase for arbitrary combinations of $`P_0`$ and $`T`$; this happens only on the cloud curve. To find the cloud curve, we combine the coexistence conditions derived above with the equation of the dilution line, so that one of the coexisting phases ($`\alpha `$ say) is the parent distribution; this implies that $`\rho ^\alpha (\sigma )=P(\sigma )`$, $`\lambda _1^\alpha =0`$, $`\rho _0^\alpha =P_0`$, and $`\rho _1^\alpha =P_0n_1`$ (or $`\rho _1^\alpha =\rho _0^\alpha n_1`$).

### D Critical Point

The critical point of the polydisperse system lies at the intersection of the cloud and shadow curves and is, therefore, correctly described (given Eq. (11)) by the projected free energy. The equations to determine it are analogous to those of the monodisperse case,

$`{\displaystyle \frac{^2F}{\rho _1^2}}|_c=0`$ (23)
and

$`{\displaystyle \frac{^3F}{\rho _1^3}}|_c=0,`$ (24)
where $`|_c`$ means evaluated at the critical point $`(\rho _1^c,T_c;P_0^c)`$. Note that these two equations have three unknowns, namely $`\rho _1^c`$, $`T_c`$ and $`P_0^c`$, and must be complemented by the dilution line constraint, $`\rho _1^c=P_0^cn_1`$.

## IV Landau Expansion

The projected free energy $`F(\rho _1,T;P_0)`$ can now be expanded about the critical point,

$`F{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2F}{\rho _1^2}}|_{\rho _1^c,T;P_0}\delta \rho ^2+`$ (25)
$`{\displaystyle \frac{1}{3!}}{\displaystyle \frac{^3F}{\rho _1^3}}|_{\rho _1^c,T;P_0}\delta \rho ^3+{\displaystyle \frac{1}{4!}}{\displaystyle \frac{^4F}{\rho _1^4}}|_{\rho _1^c,T;P_0}\delta \rho ^4+\mathrm{},`$ (26)
where $`\delta \rho \rho _1\rho _1^c`$. As before, the constant and linear terms have been dismissed.

Now we must expand the coefficients of this expansion around the critical values not only of $`T`$, as we did in the monodisperse case, but also of $`P_0`$. One might think that, after doing this, the free energy could be written, to the lowest order, as

$`F\left(A_1\delta T+A_2\delta P_0\right)\delta \rho ^2+C\delta \rho ^4,`$ (27)
where $`\delta P_0P_0P_0^c`$. If so, the term in brackets should be of order $`\delta \rho ^2`$ to balance the term of order $`\delta \rho ^4`$. However, the second term in the brackets is actually linear in $`\delta \rho `$, as dictated by the dilution line $`\delta P_0=\delta \rho /n_1`$. These two facts can only be obeyed if the first term in the brackets behaves like $`A_1\delta TA_2\delta \rho /n_1+𝒪\left(\delta \rho ^2\right)`$ (instead of $`\delta T\delta \rho ^2`$, as in the monodisperse case). This corresponds to the fact that, generically, the critical point is no longer at the top of the coexistence curve. Therefore, to maintain consistently the order of the approximation, the correct lowest order expansion of the free energy must contain all terms of order $`\delta \rho ^4`$:

$`F{\displaystyle \frac{1}{2}}[{\displaystyle \frac{^3F}{\rho _1^2T}}|_c\delta T+{\displaystyle \frac{^3F}{\rho _1^2P_0}}|_c\delta P_0+{\displaystyle \frac{1}{2}}{\displaystyle \frac{^4F}{\rho _1^2T^2}}|_c\delta T^2+`$ (28)
$`{\displaystyle \frac{1}{2}}{\displaystyle \frac{^4F}{\rho _1^2P_0^{\mathrm{\hspace{0.17em}2}}}}|_c\delta P_{0}^{}{}_{}{}^{2}+{\displaystyle \frac{^4F}{\rho _1^2TP_0}}|_c\delta T\delta P_0]\delta \rho ^2+`$ (29)
$`{\displaystyle \frac{1}{3!}}\left[{\displaystyle \frac{^4F}{\rho _1^3T}}|_c\delta T+{\displaystyle \frac{^4F}{\rho _1^3P_0}}|_c\delta P_0\right]\delta \rho ^3+`$ (30)
$`{\displaystyle \frac{1}{4!}}{\displaystyle \frac{^4F}{\rho _1^4}}|_c\delta \rho ^4+\mathrm{}.`$ (31)
To continue, let us define the functions

$`\rho _m={\displaystyle 𝑑\sigma P(\sigma )\left(w_1(\sigma )\right)^me^{\lambda _1w_1(\sigma )}},`$ (32)
for any $`m=0,1,2,\mathrm{}`$, which are an expanded set of generalised moments consistent with the previous definitions Eq. (20) and (22). These functions depend on $`\rho _1`$ and $`P_0`$ and, indeed, they obey:

$`{\displaystyle \frac{\rho _m}{\rho _1}}={\displaystyle \frac{\rho _{m+1}}{\rho _2}},`$ (33)
and

$`{\displaystyle \frac{\rho _m}{P_0}}={\displaystyle \frac{1}{P_0}}\left(\rho _m{\displaystyle \frac{\rho _1\rho _{m+1}}{\rho _2}}\right).`$ (34)
Very importantly, these equations imply that the derivatives of any (generalised) moment with respect to $`\rho _1`$ gives rise to functions of moments (Eq. (33)), and the derivatives of any moment with respect to $`P_0`$ must vanish in the monodisperse limit (as does Eq. (34)).

With these definitions, we can write the second derivative of the free energy, from which all the coefficients of the expansion (28) can be obtained, as:

$`{\displaystyle \frac{^2F}{\rho _{1}^{}{}_{}{}^{2}}}={\displaystyle \frac{1}{\rho _2}}+{\displaystyle \frac{^2f(\rho _1,T)}{\rho _{1}^{}{}_{}{}^{2}}}.`$ (35)
Now note that, since $`T`$ and $`P_0`$ appear in disjunct parts of the free energy (as in the original free energy, Eq. (21)), all crossed terms must vanish:

$`{\displaystyle \frac{^2F}{TP_0}}=0.`$ (36)
Accordingly, the coefficients of the expansion can be divided into three different groups: those including only derivatives with respect to $`\rho _1`$, those including derivatives with respect to $`T`$, and those including derivatives with respect to $`P_0`$. While the first group of coefficients originates from the full free energy, the second originates from the excess part of the free energy (‘the interactions’), and the third group arises from the ideal part of the free energy (‘the entropy’). These three groups of coefficients have a distinct dependence on polydispersity. In particular, all the coefficients stemming from the ideal part (those multiplying any power of $`\delta P_0`$) must vanish in the monodisperse limit , hence recovering the ordinary Landau expansion, Eq. (8) (plus some terms e.g. in $`\delta T^2\delta \rho ^2`$ that are higher order corrections in this limit).

In brief, the correct lowest order expansion of the free energy, retaining all terms to fourth order in density, is therefore

$`F\left(A_1\delta T+A_2\delta P_0+A_3\delta T^2+A_4\delta P_{0}^{}{}_{}{}^{2}\right)\delta \rho ^2+`$ (37)
$`\left(B_1\delta T+B_2\delta P_0\right)\delta \rho ^3+C_0\delta \rho ^4,`$ (38)
where the coefficients of terms involving $`\delta T`$ (those with odd indices: $`A_1`$, $`A_3`$ and $`B_1`$) only depend on derivatives of the excess free energy $`f`$, and the coefficients of terms involving $`\delta P_0`$ ($`A_2`$, $`A_4`$ and $`B_2`$) stem solely from the ideal entropy term of the system and not the interactions. The latter vanish in the monodisperse limit; for example,

$`A_2{\displaystyle \frac{1}{2P_{0}^{c}{}_{}{}^{\mathrm{\hspace{0.17em}2}}}}\left({\displaystyle \frac{n_1n_3}{n_{2}^{}{}_{}{}^{3}}}{\displaystyle \frac{1}{n_2}}\right),`$ (39)
where $`n_1,n_2\mathrm{}`$ are generalized moments of the (normalized) parent distribution,

$`n_m{\displaystyle 𝑑\sigma n(\sigma )\left(w_1(\sigma )\right)^m},`$ (40)
which do not depend on any thermodynamic variable.

It is interesting to point out here that the ultimate reason why the critical point is not on the top of the cloud curve is that the squared term of the free energy expansion at the critical point behaves effectively as $`\delta T^{\mathrm{\hspace{0.17em}2}}\delta \rho ^{\mathrm{\hspace{0.17em}2}}`$ instead of the usual $`\delta T\delta \rho ^{\mathrm{\hspace{0.17em}2}}`$.

We stress that although the coefficients $`A_2`$, $`A_4`$ and $`B_2`$ are small for narrow polydispersities, the expansion (37) is asymptotically exact around the critical point for arbitrary (not necessarily narrow) polydisperse distributions. As ever, ‘exact’ means exact within a mean-field free energy of the form Eq. (11).

### A Cloud Curve

In order to obtain an expression for the cloud curve, we can write the free energy expansion (37) in terms of a shifted density

$`\delta \eta \delta \rho +{\displaystyle \frac{B_1\delta T+B_2\delta P_0}{4C_0}},`$ (41)
obtaining an asymptotic expression of the form

$`FA\delta \eta ^2+C_0\delta \eta ^4,`$ (42)
where

$`AA_1\delta T+A_2\delta P_0+A_3\delta T^2+A_4\delta P_{0}^{}{}_{}{}^{2}`$ (43)
$`{\displaystyle \frac{3}{8}}{\displaystyle \frac{(B_1\delta T+B_2\delta P_0)^2}{C_0}}.`$ (44)
Starting from this simple expression, the coexisting states follow immediately as

$`\delta \eta \pm \sqrt{{\displaystyle \frac{A}{2C_0}}}.`$ (45)
These two coexisting densities correspond to distribution functions given by Eq. (19), neither of which is in general the parent distribution function $`P(\sigma )`$. In order to obtain the appropriate cloud curve, we must therefore impose the dilution line constraint, $`\delta \rho ^{CL}=n_1\delta P_0`$, which leads to the following expression for the cloud curve:

$`\delta T_{CL}\alpha _{CL}\delta \rho +\beta _{CL}\delta \rho ^2+𝒪(\delta \rho ^3),`$ (46)
where the slope of the cloud curve at the critical point is given by

$`\alpha _{CL}{\displaystyle \frac{A_2}{A_1n_1}},`$ (47)
and its curvature by

$`\beta _{CL}{\displaystyle \frac{C_0}{A_1}}\left[\left(1\mathrm{\Delta }\right)^23\right]{\displaystyle \frac{A_4}{A_1n_1^2}}{\displaystyle \frac{A_2^2A_3}{A_1^3n_1^2}},`$ (48)
where

$`\mathrm{\Delta }{\displaystyle \frac{A_1B_2A_2B_1}{2A_1C_0n_1}}.`$ (49)
As predicted (preceding Eq. (28) above), we have thus established that $`A_1\delta T+A_2\delta P_0\delta \rho ^2`$, with exact cancelation of the terms linear in $`\delta \rho `$. This shows the consistency of the approximation and proves that none of the terms in Eq. (28) can be discarded in the expansion.

One consequence of these results concerns the geometric location of the critical point along the cloud curve. For the critical point to remain at the apex of the cloud curve (which is where it resides in the monodisperse limit), the slope of the cloud curve at the critical point $`\alpha _{CL}`$, Eq. (47), must vanish; this requires $`A_2=0`$. According to Eq. (39), this requires very particular combinations of the parent distribution and weight function $`w_1(\sigma )`$. If the latter only adopts positive values, $`A_2`$ is strictly positive except in the monodisperse limit so that the critical point will move away from the top of the cloud curve even for infinitesimal polydispersity . Since the argument originates from the entropic part of the free energy, it may hold more generally than for systems obeying Eq. (11) .

For monodisperse systems, the slope $`\alpha _{CL}`$ vanishes, while $`\beta _{CL}`$ tends to $`2C_0/A_1`$, thus recovering the classical scenario that leads to Eq. (9): the cloud curve forms the binodal curve with a quadratic maximum where the critical point resides.

### B Shadow Curve

The shadow curve does not, of course, obey the dilution line constraint, and therefore is not obtainable direct from Eq. (45). Nonetheless we expect

$`\delta T_{SH}\alpha _{SH}\delta \rho +\beta _{SH}\delta \rho ^2+𝒪(\delta \rho ^3).`$ (50)
To find $`\alpha _{SH}`$, we may use Eq. (45) to write a relation between the cloud and the shadow curves,

$`\delta \eta ^{CL}+\delta \eta ^{SH}0.`$ (51)
This yields, upon substitution of Eq. (41),

$`\alpha _{SH}{\displaystyle \frac{\alpha _{CL}}{\mathrm{\hspace{0.33em}1}+\mathrm{\Delta }}}.`$ (52)
Interestingly though, Eq. (51) does not provide the correct curvature of the shadow curve, $`\beta _{SH}`$, as can easily be checked numerically for the example considered below (Section V). This is because the curvature of the shadow in fact depends on some higher-order terms (of order $`\delta \rho ^5`$) that were not included in Eq. (37). The fifth-order terms are

$`\mathrm{\Delta }F=\left(A_5\delta T^3+A_6\delta P_{0}^{}{}_{}{}^{3}\right)\delta \rho ^2+\left(B_3\delta T^2+B_4\delta P_{0}^{}{}_{}{}^{2}\right)\delta \rho ^3`$ (53)
$`+\left(C_1\delta T+C_2\delta P_0\right)\delta \rho ^4+D_0\delta \rho ^5,`$ (54)
Including these (fifth-order) terms in the expansion, Eq. (37), transforms Eq. (51) into

$`\delta \eta ^{CL}+\delta \eta ^{SH}`$ (55)
$`{\displaystyle \frac{(A_1\delta T+A_2\delta P_0+A_3\delta T^2+A_4\delta P_{0}^{}{}_{}{}^{2})D_0}{2\left(C_0+C_1\delta T+C_2\delta P_0\right)^2}}+𝒪(\delta \eta ^3),`$ (56)
from which we finally obtain the expresion

$`\beta _{SH}{\displaystyle \frac{1}{\left(\mathrm{\hspace{0.33em}1}+\mathrm{\Delta }\right)^3}}[\beta _{CL}(1+{\displaystyle \frac{1}{2C_0^2n_1}}(B_2C_0A_2D_0))`$ (57)
$`+{\displaystyle \frac{\alpha _{CL}}{2C_0^2n_1^2}}(B_2C_2B_4C_0+A_4D_0)`$ (58)
$`+{\displaystyle \frac{\alpha _{CL}^{\mathrm{\hspace{0.17em}2}}}{2C_0^2n_1}}(B_1C_2+B_2C_1)`$ (59)
$`+{\displaystyle \frac{\alpha _{CL}^{\mathrm{\hspace{0.17em}3}}}{2C_0^2}}(B_1C_1B_3C_0+A_3D_0)].`$ (60)
In the monodisperse limit, analogously to the cloud curve, the slope $`\alpha _{SH}`$ vanishes and $`\beta _{SH}`$ tends to $`2C_0/A_1`$, thus forming the binodal curve by the merger of cloud and shadow.

### C Consequences of the Expansion

Note the progression in complexity of the coefficients needed in the expansion to compute the slopes, and then curvatures, of the cloud and shadow at the critical point. The slope of the cloud curve $`\alpha _{CL}`$ depends on just two coefficients of the expansion ($`A_1`$ and $`A_2`$); the slope of the shadow involves three extra coefficients ($`B_1`$, $`B_2`$ and $`C_0`$); the curvature of the cloud at the critical point requires knowledge of all seven coefficients of the quartic Landau expansion, Eq. (37); that of the shadow also involves a selection of the fifth order terms (but not all of them): Table 1 shows which terms in the Landau expansion must be retained to correctly predict each of these parameters.

Although these results are obtained for a one-moment free-energy, Eq. (11), they involve higher moments of the parent distribution , as detailed in Table 1. Note also that the critical point and the spinodals depend only on a reduced set of moments; this holds to all orders in the Landau expansion. In contrast, order by order, the cloud and shadow curves involve progressively all the moments of the parent distribution, through Eqs. (33) and (34).

The results given above were obtained by the projection method, but since this gives the exact cloud and shadow curves, any other method for finding $`\alpha _{CL}`$, $`\beta _{CL}\mathrm{}`$ within mean field theory must produce, necessarily, expressions equivalent to Eqs. (47,48,52,57).

### D Phase Diagram in Temperature-Pressure Plane

In the monodisperse case, the phase equilibria show up in the $`Tp`$ phase diagram as a line of coexistence, ending at the critical point. If polydispersity is included, this line now broadens and adopts a banana-like shape of finite area, containing the values of temperature and pressure for which the system undergoes phase separation. In this case, the critical point appears on the perimeter of the coexistence region.

We now can calculate the shape of the coexistence region to the same order of approximation as we did with the cloud and the shadow. After projection of the free energy to give Eq. (21), the pressure of the system can be written in the familiar form

$`{\displaystyle \frac{p}{k_BT}}={\displaystyle \frac{F}{\rho _1}}\rho _1F,`$ (61)
which allows us to expand $`\delta ppp_c`$ as

$`{\displaystyle \frac{\delta p}{k_BT}}a_1\delta T+a_2\delta P_0+b_1\delta T^2+b_2\delta P_{0}^{}{}_{}{}^{2}+\mathrm{}`$ (62)
where (just as occurred in Eq. (37)) the coefficients of terms involving $`\delta T`$ (those with odd indices: $`a_1`$ and $`b_1`$) only depend on derivatives of the excess free energy $`f`$, and the coefficients of terms involving $`\delta P_0`$ ($`a_2`$ and $`b_2`$) derive solely from the ideal entropy of the system (not the interactions) and vanish in the monodisperse limit. For instance,

$`a_2{\displaystyle \frac{n_1^{\mathrm{\hspace{0.17em}2}}}{n_2}}1(0).`$ (63)
Substitution of the approximation derived in Section IV A for the cloud curve yields

$`{\displaystyle \frac{\delta p}{k_BT}}\alpha _P\delta \rho +\beta _P\delta \rho ^2,`$ (64)
where

$`\alpha _Pa_1\alpha _{CL}+{\displaystyle \frac{a_2}{n_1}}`$ (65)
and

$`\beta _Pa_1\beta _{CL}+b_1\alpha _{CL}^{}{}_{}{}^{2}+{\displaystyle \frac{b_2}{n_1^2}}.`$ (66)
A generic consequence of these results concerns the location of the critical point along the perimeter of the coexistence region on the $`TP`$ phase diagram. It is clear on geometrical grounds that the critical point cannot be located at the point of maximum temperature and maximum pressure simultaneously, unless the coexistence region shows a kink at that point. Such a kink is not predicted by the Landau expansion, which gives smooth behavior everywhere; nor have we seen evidence for it in any of the models (or experiments) studied to date. Assuming these two turning points are indeed distinct, one can ask whether the critical point generically lies at one, at the other, or neither. We argued already above that generically, for systems obeying Eq. (11), the critical point is not at the maximum of $`T`$ since this would require $`A_2`$ to vanish (Eq. 39). Likewise, for it to lie at the pressure maximum would require $`\alpha _P=0`$ in Eq. (65), for which a fine tuning of the parameters appears necessary . We conclude that generically it lies at neither of these points; indeed, as shown in the following example, it lies between them.

## V Example

In this Section, we illustrate the results of the Landau expansion for the cloud and shadow curves of a particular system: an incompressible blend of two homopolymers. One of the polymers is considered polydisperse in length, with number densities $`\rho (\sigma )`$ where $`\sigma `$ is the number of monomers on a chain (effectively a continuous parameter). The other polymer (or solvent) is monodisperse and has $`N^{}`$ monomers. For this system, the free energy can be easily obtained from the Flory-Huggins theory and recovers the form (11). Specifically, the excess free energy reads

$`f(\rho _1,\chi )={\displaystyle \frac{1\rho _1}{N^{}}}\mathrm{log}\left({\displaystyle \frac{1\rho _1}{N^{}}}\right)+\chi \rho _1(1\rho _1),`$ (67)
where $`\chi `$ is the Flory parameter, which plays the role of an inverse temperature. The weight function for the moment $`\rho _1`$ is $`w_1(\sigma )=\sigma `$ (in the usual polymer jargon, $`\rho _1`$ is the monomer volume fraction $`\varphi `$).

As a normalised parent, we consider in this example a uniform distribution

$$n(\sigma )=\{\begin{array}{cc}\frac{1}{\mathrm{\hspace{0.33em}2}\sqrt{3}\mathrm{\Delta }N}& \text{ if }|\sigma N|\sqrt{3}\mathrm{\Delta }N\hfill \\ & \\ 0& \text{ otherwise}\hfill \end{array}$$
(68)
where $`N=n_1=\sigma `$ is the mean polymer size, and $`\mathrm{\Delta }N=(\sigma ^2\sigma ^2)^{1/2}`$ is the standard deviation.

We plot in Figures 1 and 2 the resulting behavior for the case where the polymer solvent has $`N^{}=10`$ monomers and the distribution of sizes of the solute follows Eq. (68) with $`N=1000`$ and $`\mathrm{\Delta }N/N=0.2`$. (This corresponds, in polymer jargon, to a polydispersity index $`M_w/M_n=\sigma ^2/\sigma ^2=1.4`$.)

We see from these plots that the Landau expansion, Eq. (37), captures correctly all the essential features of the full mean-field model near the critical region, just as it does for monodisperse systems. In this example, the Landau predictions cross the cloud and shadow curves at the critical point (although the slopes and curvatures are correct). This behavior is inherited from the monodisperse case, and is presumably connected with an error in the third derivative of these curves.

Given the increased complexity of the phase behavior (separation of cloud and shadow curves; critical point not at the maximum of these; finite width coexistence region on the $`pT`$ plot collapsing to a line in the monodisperse limit), we find this somewhat remarkable. We emphasize that every term in Eq. (37) must be retained to achieve this degree of consistency – and in particular to reproduce the correct curvature of the cloud; to get the curvature of the shadow also correct, one needs as well the higher terms in Eq. (53). The latter were included in the figure although leaving them out gives a result that, for this example, is only slightly less accurate. Despite all the extra complexity, the region of the phase diagram within which the Landau expansion could tolerably be used to replace the full calculation is similar to the monodisperse case.

## VI Conclusions

Within the one-moment form of the excess free energy represented by Eq. (11), we have derived a generalized Landau expansion, Eqs. (37,53), from which the characteristic features of polydisperse binary fluid coexistence can all be shown to follow. This Landau expansion was compared with direct calculation from the underlying (mean-field) free energy in a specific example (Figures 1 and 2) with excellent agreement near the critical point and qualitative agreement away from it. Note that truncation of the Landau expansion at the quartic order (Eq.37) may prove adequate for qualitative purposes if the curvature of the shadow is not important.

Our expansion is correct around the critical point for arbitrary polydispersity, be it broad or narrow. The monodisperse limit is recovered by decreasing the polydispersity. In a pressure-temperature representation this limit is highly singular: for example the slope of the perimeter of the coexistence region at the critical point will converge, as the monodisperse limit is taken, to a negative value – even though the coexistence line (onto which the region itself collapses in this limit) has positive slope everywhere! This behavior is easily understood from Figure 2, and the Landau expansion seems to have no difficulty dealing with it.

The perturbation series in density that underlies the Landau expansion contrasts with, and is complementary to, a previous perturbative approach that was developed for the limit of narrow polydisperse distributions . The latter approach always breaks down close enough to the critical point (even within mean field theory, to which, unlike our work, it is not limited in principle).

Finally, we note that although our Landau expansion was derived within the specific context of Eq. (11), the generic features of the phase diagram that it reproduces are seen rather generally in experiments and other theoretical models. We have not attempted explicit calculation for the case where more than one moment is present in the excess free energy, but since the topology of the phase diagram (near the fluid-fluid critical point) is already captured by Eq. (11), it is possible that any such refinements will not alter the basic structure of these results. A firm conclusion on this point is, however, left to future work.

Acknowledgements We are grateful to Yuri Martínez-Ratón, José Cuesta, Patrick Warren and Peter Sollich for discussions. Work funded by EPSRC Grant GR/M29696. C.R. also acknowledges support by the Ministerio de Ciencia y Tecnología (Spain) under grant BFM2000-0004.


text:
# IFIC/02-43, FTUV/02-0930 NNLO 𝜏⁺𝜏⁻ production cross section at threshold

## 1 Introduction

The Tau–Charm Factory, a high–luminosity ($`10^{33}\text{cm}^2\text{s}^1`$) $`e^+e^{}`$ collider with a centre–of–mass energy near the $`\tau ^+\tau ^{}`$ production threshold, has been proposed as a powerful tool to perform high–precision studies of the $`\tau `$ lepton, charm hadrons and the charmonium system . In recent years, this energy region has been only partially explored by the Chinese BEBC machine ($`10^{31}\text{cm}^2\text{s}^1`$). The possibility to operate the Cornell CESR collider around the $`\tau ^+\tau ^{}`$ threshold has revived again the interest on Tau–Charm Factory physics.

A precise understanding of the $`e^+e^{}\tau ^+\tau ^{}`$ production cross section near threshold is clearly required. The accurate experimental analysis of this observable could allow to improve the present measurement  of the $`\tau `$ lepton mass. The cross section $`\sigma (e^+e^{}\tau ^+\tau ^{})`$ has already been analysed to $`𝒪(\alpha ^3)`$ in refs. , including a resummation of the leading Coulomb corrections.

The recent development of non-relativistic effective field theories of QED (NRQED) and QCD (NRQCD)  has allowed an extensive investigation of the threshold production of heavy flavours at $`e^+e^{}`$ colliders. The threshold $`b\overline{b}`$  and $`t\overline{t}`$  production cross sections have been computed to the next-to-next-to-leading order (NNLO) in a combined expansion in powers of $`\alpha _s`$ and the fermion velocities. Making appropriate changes, those calculations can be easily applied to the study of $`\tau ^+\tau ^{}`$ production . One can then achieve a theoretical precision better than 0.1%.

## 2 Perturbative calculation to $`𝒪(\alpha ^4)`$

At lowest order in QED, the $`\tau `$ leptons are produced by one-photon exchange in the s-channel, and the total cross section formula reads

$$\sigma _B(e^+e^{}\tau ^+\tau ^{})=\frac{2\pi \alpha ^2}{3s}v(3v^2),$$
(1)
where $`v=\sqrt{14m_\tau ^2/s}`$ is the velocity of the final $`\tau `$ leptons in the center-of-mass frame of the $`e^+e^{}`$ pair which makes $`\sigma _B`$ vanish when $`v0`$.

Electromagnetic corrections of $`𝒪(\alpha )`$ arise from the interference between the tree level result and 1-loop amplitudes. A factor $`\alpha /v`$ emerges in the 1-loop final state interaction between the tau leptons, making the cross section at threshold finite. Furry’s theorem guarantees that contributions to $`\sigma (e^+e^{}\tau ^+\tau ^{})`$ coming from initial, intermediate and final state corrections completely factorize at $`𝒪(\alpha ^3)`$, including real photon emission.

Some undesirable features appear at $`𝒪(\alpha ^4)`$: The two-loop $`\tau ^+\tau ^{}\gamma `$ vertex develops an $`\alpha ^2/v^2`$ term which makes the cross section ill-defined when $`v0`$, and multiple photon production of tau leptons by box-type diagrams and the non-zero interference of initial and final state radition spoil exact factorization. However, as it has been recently shown , the squared amplitude of the $`e^+e^{}\tau ^+\tau ^{}`$ box diagram is proportional to $`\alpha ^4v^2`$, and so represents a N4LO correction in the combined expansion in powers of $`\alpha `$ and $`v`$, far beyond the scope of this analysis. In addition, contributions to the total cross section from diagrams with real photons emitted from the produced taus can be shown to begin at N3LO, and factorization remains at NNLO. The total cross section can thus be written as an integration over the product of separate pieces including initial, intermediate and final state corrections:

$$\sigma (s)=^sF(s,w)\left|\frac{1}{1+e^2\mathrm{\Pi }\text{em}(w)}\right|^2\stackrel{~}{\sigma }(w)𝑑w.$$
(2)
The radiation function $`F(s,w)`$  describes initial state radiation, including virtual corrections. The integration accounts for the effective energy loss due to photon emission from the $`e^+e^{}`$ pair, and it includes the largest corrections coming from the emission of an arbitrary number of initial photons, which can sizeably suppress the total cross section. $`\stackrel{~}{\sigma }`$ collects only final-state interactions between the tau leptons, and it is usually written in terms of the tau spectral density $`R_\tau `$,

$$\stackrel{~}{\sigma }(e^+e^{}\gamma ^{}\tau ^+\tau ^{})=R_\tau (s)\sigma _{pt},$$
(3)
with $`\sigma _{pt}=\frac{4\pi \alpha ^2}{3s}`$. The threshold behaviour of the total cross section will be ruled by the expansion of $`R_\tau `$ at low velocities.

## 3 Non-Relativistic Corrections: NRQED

A NNLO calculation of the cross section in the kinematic region where $`v\alpha `$ has to account for all terms proportional to $`v(\alpha /v)^n\times [1;\alpha ;v;\alpha ^2;\alpha v;v^2]`$ with $`n=1,2,\mathrm{}`$ The leading divergences (i.e. $`\left(\frac{\alpha }{v}\right)^n`$, $`n>1`$) are resummed in the well-known Sommerfeld factor

$$|\mathrm{\Psi }_{c,\text{E}}(0)|^2=\frac{\alpha \pi /v}{1\text{exp}(\alpha \pi /v)},$$
(4)
multiplying the Born cross section (1). A systematic way to calculate higher-order corrections in this regime requires the use of a simplified theory which keeps the relevant physics at the scale $`MvM\alpha `$, characteristic of the Coulomb interaction. NRQED  was designed precisely for this purpose. It is an effective field theory of QED at low energies, applicable to fermions in non-relativistic regimes, i.e. with typical momenta $`p/Mv1`$. Interactions contained in the NRQED Lagrangian have a definite velocity counting but propagators and loop integrations can also generate powers of $`v`$. With appropriate counting rules at hand, one can prove that all interactions between the non-relativistic pair $`\tau ^+\tau ^{}`$ can be described up to NNLO in terms of time-independent potentials , derived from the low-energy Lagrangian. Therefore, the low-energy expression of the $`\tau `$ spectral density is related with the non-relativistic Green’s functions :

$$R_\tau ^{\text{NNLO}}(s)=\frac{6\pi }{M^2}\text{Im}\left(C_1G(E)\frac{4E}{3M}G_c(E)\right)$$
(5)
with $`C_1`$ a short distance coefficient to be determined by matching full and effective theory results and $`E=m_\tau v^2`$ the non-relativistic energy. The details of this derivation can be found in the Appendix B of . The Green’s function $`G`$ obeys the Schrödinger equation corresponding to a two-body system interacting through potentials derived from $`_{\text{NRQED}}`$ at NNLO:

$`({\displaystyle \frac{\mathbf{}^2}{M}}`$ $``$ $`{\displaystyle \frac{\mathbf{}^4}{4M^3}}+V_c(𝒓)+V_{\text{BF}}(𝒓)+V_{\text{An}}(𝒓)E)`$ (6)
$`\times `$ $`G(𝒓,𝒓^{},E)=\delta ^{(3)}(𝒓𝒓^{})`$
The term $`\frac{\mathbf{}^4}{4M^3}`$ is the first relativistic correction to the kinetic energy. $`V_c`$ stands for the Coulomb potential with $`𝒪(\alpha ^2)`$ corrections. At NNLO, the heavy leptons are only produced in triplet S-wave states, so we just need to consider the corresponding projection of the Breit-Fermi potential $`V_{\text{BF}}`$. Finally, $`V_{\text{An}}`$ is a NNLO piece derived from a contact term in $`_{\text{NRQED}}`$, which reproduces the QED tree level s-channel diagram for the process $`\tau ^+\tau ^{}\tau ^+\tau ^{}`$.

A solution of eq. (6) must rely on numerical or perturbative techniques. In the QED case, a significant difference between both approaches is not expected, being $`\alpha `$ such a small parameter. Consequently we follow the perturbative approach, using recent results , where NLO and NNLO corrections to the Green’s function are calculated analytically, via Rayleigh-Schrödinger time-independent perturbation theory around the known LO Coulomb Green’s function $`G_c`$. We refer the reader to Appendix C of  for complete expressions of the Green’s function corrections.

## 4 Vacuum Polarization

For a complete NNLO description of $`\sigma (e^+e^{}\tau ^+\tau ^{})`$, two-loop corrections to the photon propagator should be included. The light lepton contributions to the vacuum polarization are the standard 1- and 2-loop perturbative expressions. For the $`\tau `$ contribution in the threshold vicinity $`q^2\stackrel{>}{_{}}4M^2`$, resummation of singular terms in the limit $`v0`$ is again mandatory. Under the assumption $`\alpha v`$, we need to know NLO contributions to $`\mathrm{\Pi }_\tau (q^2)`$, performing the direct matching for both real and imaginary parts.

In the hadronic sector, we can relate the hadronic vacuum polarization with the total cross section $`\sigma (e^+e^{}\gamma ^{}had)`$. Below 1 GeV, the electromagnetic production of hadrons is dominated by the $`\rho `$ resonance and its decay to two charged pions. The photon mediated $`\pi ^+\pi ^{}`$ production cross section is driven by the pion electromagnetic form factor $`F(s)`$. An analytic expression for $`F(s)`$ was obtained in Ref. , using Resonance Chiral Theory and the restrictions imposed by analyticity and unitarity. The obtained $`F(s)`$ provides an excellent description of experimental data up to energies of the order $`s_\rho `$ 1 GeV2.

For the integration region above $`s_\rho `$, we use $`\text{Im}\mathrm{\Pi }_{\text{had}}`$ as calculated from pQCD. Our simple estimate has been proved  to deviate by less than 5$`\%`$ for the running of $`\alpha `$ at the scale $`\sqrt{s}=M_Z`$. Considering that $`\mathrm{\Pi }_{\text{had}}`$ modifies $`\sigma (e^+e^{}\tau ^+\tau ^{})`$ near threshold by roughly $`1\%`$, our result has a global uncertainty smaller than $`0.1\%`$ for the total cross section.

## 5 Bound states effect

Green’s functions develop energy poles below threshold corresponding to spin triplet (n$`{}_{}{}^{3}S_{1}^{}`$) electromagnetic $`\tau ^+\tau ^{}`$ bound states. The small width of these bound levels, dominated by their $`e^+e^{}`$ decay rate

$$\mathrm{\Gamma }_{ee}=\frac{m_\tau \alpha ^5}{6n^3}\frac{6.110^3\text{eV}}{n^3},$$

make these states very difficult to be resolved experimentally. For the same reason, the bound states cannot affect the shape of the cross section at threshold, contrary to the case of heavy quark threshold production, where bound states play a crucial role.

## 6 Numerical analysis for $`\sigma (e^+e^{}\tau ^+\tau ^{})`$

The need for performing resummations of the leading non-relativistic terms $`\left(\alpha /v\right)^n[v,v\alpha ,v^2,\mathrm{}]`$ is evidenced in Fig. 1. The QED spectral density vanishes as $`v0`$, due to the phase space velocity in formula (1), which is cancelled by the first $`v^1`$ term appearing in the $`𝒪(\alpha )`$ correction, making the cross section at threshold finite. More singular terms near threshold, $`v^2,\mathrm{}`$ arising in higher-order corrections completely spoil the expected good convergence of the QED perturbative series at low $`v`$. This is no longer the case for the effective theory perturbative series, whose convergence improves as we approach the threshold point, as shown in Fig. 2, and higher-order corrections reduce the perturbative uncertainty. In the whole energy range displayed in Fig. 2, the differences between the NNLO, NLO and LO results are below 0.8%, which indicates that the LO result, i.e. the Sommerfeld factor, contains the relevant physics to describe the threshold region.

Adding the intermediate and initial state corrections we have a complete description of the total cross section of $`\tau ^+\tau ^{}`$ production, as shown in Fig. 3. Coulomb interaction between the produced $`\tau `$’s, becomes essential within few MeV’s above the threshold, and its effects have to be taken into account to all orders. Initial state radiation effectively reduces the available center-of-mass energy for $`\tau `$ production, lowering in this way the total cross section.

We should emphasize that NNLO corrections do not modify the predicted behaviour of the LO and NLO cross section, but are essential to guarantee that the truncated perturbative series at NLO gets small corrections from higher-order terms. Hence, we have shown that the theoretical uncertainty of our analysis of $`\sigma (e^+e^{}\tau ^+\tau ^{})`$ is lower than 0.1%. Nevertheless, the statistical uncertainty of the most recent experiments is still much larger than the theoretical one due to low statistics, and we should wait for future machines to improve it.

Acknowledgements I wish to thank S. Narison for the organization of the QCD 02 conference. This work has been supported in part by TMR, EC Contract No. ERB FMRX-CT98-0169, by MCYT (Spain) under grant FPA2001-3031, and by ERDF funds from the European Commission.