format: text

text:
# Spin foam quantization of SO(4) Plebanski’s action

## 1 Introduction

In reference  Barrett and Crane introduced a very interesting model of Riemannian quantum gravity based on a constrained state-sum. The definition of the model can be nicely motivated by geometrical properties of the so-called ‘quantum tetrahedron’. The definition of the quantum tetrahedron in 3 dimensions was originally introduced by Barbieri in . Baez and Barrett showed that a generalization to 4 dimensions naturally leads to the the Barrett-Crane (BC) model . Evidence suggesting that the model corresponds to a discrete path integral for general relativity has been found in . The model turns out to be well defined on a finite (non-degenerate) triangulation once an appropriate normalization is chosen. This normalization arises naturally in the so-called group field theory (GFT) formulation , which in addition provides a prescription for summing over discretizations. The model has been extended to the Lorentzian sector in . The finiteness properties are preserved in this extension.

The $`SO(4)`$ Plebanski action corresponds to the $`SO(4)`$ BF action plus certain Lagrange multiplier terms imposing constraints on the $`B`$ field. Therefore, one can formally quantize the theory restricting the BF-path-integral to paths that satisfy the B-field constraints. In the literature, there is an implicit assumption that the BC model corresponds to a realization of this idea. In other words, the definition of the quantum tetrahedron in 4d (giving rise to the BC model) is sometimes regarded as an alternative way to impose the required restrictions on the B-configurations in the discretization. The purpose of this work is to analyze if this is the case by systematically carrying out this restrictions.

We will present a construction which defines the path integral of Plebanski’s action on a fixed simplicial decomposition of space-time. As just mentioned, this is done by appropriately restricting the state-sum of the $`SO(4)`$ BF theory. The path-integral of the BF theory is defined on a triangulation using techniques similar to those in lattice-gauge theory. The spin-foam formulation —or state-sum— is obtained by performing the mode expansion of certain distributions on $`SO(4)`$. This is analogous to a Fourier transform where modes correspond to unitary irreducible representations of $`SO(4)`$ (Peter-Weyl theorem). The constraints on the B-field in the classical action can be naturally translated into restrictions on these modes. The definition of these constraints is not different in spirit from that of Barrett and Crane. However, we emphasize the requirement that the restrictions should be imposed on configurations of the BF theory. After making some natural definitions, a systematic derivation leads to a model that is closely related to the BC model but that does not agree with it. This new version has a puzzling feature: states of 3-geometries (boundary spin-network states) are annihilated by the volume operator. The point of view is related to that of Reisenberger and Freidel-Krasnov in . No obvious modification of the prescription can lead to all the BC configurations.

In order to find a possible interpretation of this result we concentrate on one of the degenerate sectors of Plebanski’s action described in . It turns out that one can define a spin foam quantum model corresponding to this sector in a straightforward way. For this, one simply applies the same techniques used in the case of the BF theory. Surprisingly, the model obtained coincides with the one introduced by De Pietri et al. in . This model was defined as an alternative to the BC model arising naturally in the context of the group field theory (GFT) framework. Our result provides a clear-cut interpretations of the De Pietri et al. formulation as a quantization of a classical action. An interesting result is that all the allowed configurations for a 4-simplex in the previous model (corresponding to generic theory) are special configurations of this model.

The article is organized in the following way. In the next section we recall essential facts about $`SO(4)`$ Plebanski formulation. In Section 3 we briefly review the spin foam quantization of the BF theory and introduce our basic definitions. In Section 4 we solve the constraints that lead one from the BF theory to general relativity and construct the corresponding state-sum model. We interpret the results and show that configurations have zero 3-volume. In Section 5 we quantize the effective action corresponding to the degenerate sectors of Plebanski’s action and show that the previous model corresponds to a sub-set of the spin foam configurations obtained in the degenerate sector. We end with concluding remarks in Section 6.

## 2 Classical $`SO(4)`$ Plebanski action

Let us start by briefly reviewing Plebanski’s formulation at the classical level. Plebanski’s Riemannian action depends on an $`so(4)`$ connection $`A`$, a Lie-algebra-valued 2-form $`B`$ and Lagrange multiplier fields $`\lambda `$ and $`\mu `$. Writing explicitly the Lie-algebra indices, the action is given by

$$S[B,A,\lambda ,\mu ]=\left[B^{IJ}F_{IJ}(A)+\lambda _{IJKL}B^{IJ}B^{KL}+\mu ϵ^{IJKL}\lambda _{IJKL}\right],$$
(1)
where $`\mu `$ is a 4-form and $`\lambda _{IJKL}=\lambda _{JLKI}=\lambda _{IJLK}=\lambda _{KLIJ}`$ is tensor in the internal space. Variation with respect to $`\mu `$ imposes the constraint $`ϵ^{IJKL}\lambda _{IJKL}=0`$ on $`\lambda _{IJKL}`$. $`\lambda _{IJKL}`$ has then $`20`$ independent components. Variation with respect to the Lagrange multiplier $`\lambda `$ imposes $`20`$ algebraic equations on the $`36`$ $`B`$. Solving for $`\mu `$ they are

$$B^{IJ}B^{KL}\frac{1}{4!}ϵ_{OPQR}B^{OP}B^{QR}ϵ_{IJKL}=0$$
(2)
which is equivalent to

$$ϵ_{IJKL}B_{\mu \nu }^{IJ}B_{\rho \sigma }^{KL}=eϵ_{\mu \nu \rho \sigma },$$
(3)
for $`e0`$ where $`e=\frac{1}{4!}ϵ_{OPQR}B_{\mu \nu }^{OP}B_{\rho \sigma }^{QR}ϵ^{\mu \nu \rho \sigma }`$. The solutions to these equations are

$$B=\pm {}_{}{}^{}(ee),\mathrm{and}B=\pm ee,$$
(4)
in terms of the $`16`$ remaining degrees of freedom of the tetrad field $`e_a^I`$. If one substitutes the first solution into the original action one obtains an effective action that is precisely that of general relativity in the Palatini formulation

$$S[e,A]=\mathrm{Tr}\left[ee{}_{}{}^{}F(A)\right].$$
(5)
## 3 Quantum $`SO(4)`$ BF theory

Classical ($`Spin(4)`$) BF theory is defined by the action

$$S[B,A]=\mathrm{Tr}\left[BF(A)\right],$$
(6)
where $`B_{ab}^{IJ}`$ is a $`Spin(4)`$ Lie-algebra valued 2-form, $`A_a^{IJ}`$ is a connection on a $`Spin(4)`$ principal bundle over $``$. The theory is rather trivial and all classical solutions are locally equivalent (up to gauge transformations). The theory has only global degrees of freedom.

One can quantize the theory à la Feynman introducing a path integral measure. This is easily done by replacing the manifold $``$ by an arbitrary simplicial decomposition $`\mathrm{\Delta }`$ 111More generally, the path integral for the BF theory can be defined on an arbitrary cellular decomposition of $``$. See .. Take a fixed triangulation $`\mathrm{\Delta }`$ of $``$. The 2-skeleton of the dual of the triangulation defines a cellular 2-complex $`\mathrm{\Delta }^{}`$. Associate $`B_fso(4)`$ to each triangle in $`\mathrm{\Delta }`$ (for convenience we use the face sub index $`f`$ since triangles are in one-to-one correspondence to faces $`f\mathrm{\Delta }^{}`$), and a group element $`g_eSpin(4)`$ to each edge $`e\mathrm{\Delta }^{}`$. Consider the holonomy around faces $`U_f=g_{e_1}g_{e_2}\mathrm{}g_{e_n}`$, i.e., the product of group elements of the corresponding edges around one face (an arbitrary orientation of faces has been chosen). The discretized version of the partition function becomes

$$𝒵(\mathrm{\Delta })=\underset{f\mathrm{\Delta }^{}}{}dB_f^{\left(6\right)}\underset{e\mathrm{\Delta }^{}}{}dg_ee^{i\mathrm{Tr}\left[B_fU_f\right]}.$$
(7)
The measure $`dB_f^{\left(6\right)}`$ is the Lebesgue measure on $`\mathrm{}^6`$, while $`dg`$ corresponds to the normalized Haar measure of $`Spin(4)`$. Now the integration over the $`B_f`$’s can be done explicitly , and the result is:

$$Z(\mathrm{\Delta })=\underset{e\mathrm{\Delta }^{}}{}dg_e\underset{f\mathrm{\Delta }^{}}{}\delta (g_{e_1}\mathrm{}g_{e_n}).$$
(8)
Expanding the delta distribution in unitary irreducible representations (Peter-Weyl decomposition 222 Peter-Weyl theorem implies that
$$\delta (g)=\underset{\rho }{}\mathrm{\Delta }_\rho \mathrm{Tr}[\rho (g)],$$
where $`\rho (g)`$ is the unitary irreducible representation of dimension $`\mathrm{\Delta }_\rho `$.) we obtain

$$𝒵(\mathrm{\Delta })=\underset{𝒞:\{\rho \}\{f\}}{}\underset{e\mathrm{\Delta }^{}}{}dg_e\underset{f\mathrm{\Delta }^{}}{}\mathrm{\Delta }_{\rho _f}\mathrm{Tr}\left[\rho _f(g_e^1\mathrm{}g_e^N)\right],$$
(9)
where $`𝒞:\{\rho \}\{f\}`$ denotes the assignment of irreducible representations to faces in the dual 2-complex $`\mathrm{\Delta }^{}`$. Each particular assignment is referred to as a coloring, $`𝒞`$.

Next step is to integrate over the connection $`g_e`$. Since edges $`e\mathrm{\Delta }^{}`$ bound four different faces, each group element $`g_e`$ appears in the mode expansion of four delta functions in (8). The formula we need is that of the projection operator into the trivial component of the tensor product of four irreducible representations, namely

$$𝑑g\rho _1(g)\rho _2(g)\rho _3(g)\rho _3(g)=\underset{\iota }{}C_{\rho _1\rho _2\rho _3\rho _4}^\iota C_{\rho _1\rho _2\rho _3\rho _4}^\iota ,$$
(10)
where $`C^\iota _{\rho _1}_{\rho _2}_{\rho _3}_{\rho _4}`$ represents an orthonormal basis of invariant vectors and the sum on the RHS ranges over all the basis elements $`\iota `$.

The RHS of equation (12) can be represented graphically as

$`{\displaystyle \underset{\iota }{}}C_{\rho _1\rho _2\rho _3\rho _4}^\iota C_{\rho _1\rho _2\rho _3\rho _4}^\iota ={\displaystyle \underset{\iota }{}}\begin{array}{c}\text{}\end{array}.`$ (12)
Therefore, integrating over the $`g_e`$’s using (12) and keeping track of indices we obtain

$`Z_{BF}(\mathrm{\Delta })={\displaystyle \underset{𝒞_f:\{f\}\rho _f}{}}{\displaystyle \underset{𝒞_e:\{e\}\{\iota _e\}}{}}{\displaystyle \underset{f\mathrm{\Delta }^{}}{}}\mathrm{\Delta }_\rho {\displaystyle \underset{v\mathrm{\Delta }^{}}{}}\begin{array}{c}\text{}\end{array},`$ (14)
where the pentagonal diagram representing the vertex amplitude denotes the trace of the product of five intertwiners $`C_{\rho _1\rho _2\rho _3\rho _4}^\iota `$ according to the graphical notation of (12). Vertices $`v\mathrm{\Delta }^{}`$ are in one-to-one correspondence to 4-simplexes in the triangulation $`\mathrm{\Delta }`$. In addition we also have $`𝒞_e:\{e\}\{\iota _e\}`$ representing the assignment of intertwiners to edges. The sum over the coloring of edges, $`𝒞_e`$, comes from (12) (for an extensive explanation of the construction of the state-sum for the BF theory and the notation used here see).

What happened in going from equation (7) to (9)? We have replaced the continuous multiple integral over the $`B`$’s by the sum over representations of $`SO(4)`$. Roughly speaking, the degrees of freedom of $`B`$ are now encoded in the representation being summed over in (9). One can make a more precise definition of what ‘$`B`$’ is at the level of (9). In order to motivate our definition we isolate a single face contribution to the integrand in the partition function (7). Then we notice that the right invariant vector field $`i𝒳^{IJ}(U)`$ has a well defined action at the level of equation (9) and acts as a ‘quantum’ B at the level of (7) since

$`i𝒳^{IJ}(U)\left(e^{i\mathrm{Tr}[BU]}\right)|_{U1}=X_\nu ^{IJ\mu }U_\sigma ^\nu {\displaystyle \frac{}{U_\sigma ^\mu }}e^{i\mathrm{Tr}[BU]}|_{U1}=`$ (15)
$`=\mathrm{Tr}[X^{IJ}UB]e^{i\mathrm{Tr}[BU]}|_{U1}B^{IJ}e^{i\mathrm{Tr}[BU]},`$
where $`X^{IJ}`$ are elements of an orthonormal basis in the $`SO(4)`$ Lie-algebra. The evaluation at $`U=1`$ is motivated by the fact that configurations in the BF partition function (8) have support on flat connections.

The constraints in (3) are quadratic in the $`B`$’s. We have then to worry about cross terms, more precisely the nontrivial case corresponds to:

$`ϵ_{IJKL}𝒳^{IJ}(U)𝒳^{KL}(U)\left(e^{i\mathrm{Tr}[BU]}\right)|_{U1}`$
$`=ϵ_{IJKL}\left(\mathrm{Tr}[X^{IJ}UB]\mathrm{Tr}[X^{KL}UB]e^{i\mathrm{Tr}[BU]}+i\mathrm{Tr}[X^{IJ}X^{KL}UB]e^{i\mathrm{Tr}[BU]}\right)|_{U1}`$
$`ϵ_{IJKL}B^{IJ}B^{KL}e^{i\mathrm{Tr}[BU]},`$ (16)
where the second term on the second line can be dropped using that $`ϵ_{IJKL}X^{IJ}`$
$`X^{KL}1`$ (one of the two $`SO(4)`$ Casimir operators) and $`U1`$. Therefore, we define the $`B_f`$ field associated to a face at the level of equation (9) as the appropriate right invariant vector field $`i𝒳^{IJ}(U_f)`$ acting on the corresponding discrete holonomy $`U_f`$, namely

$$B_f^{IJ}i𝒳^{IJ}(U_f).$$
(17)
It is easy to verify that one can use left invariant vector fields instead in the previous definition without changing the following results.

## 4 Implementation of the constraints that reduce the BF theory to general relativity

### 4.1 Formulation of the problem

Now we describe the implementation of the constraints (3). The idea is to concentrate on a single 4-simplex amplitude using the locality of the BF theory state sum 333The term ‘local’ here is used as defined by Reisenberger in . It means that the spin foam can be written as 4-simplex contributions that communicate with other 4-simplexes by boundary data (connection). The full amplitude is obtained by integrating out the boundary connections along the common boundary of the 4-simplexes that make up the simplicial complex.. The 4-simplex wave function is obtained using (8) on the dual 2-complex with boundary defined by the intersection of the dual of a single 4-simplex with a 3-sphere, see Figure 1. We refer to this fundamental building block as ‘atom’ as in .

The boundary values of the discrete connection are held fixed. We denote as $`h_{ij}Spin(4)`$ ($`ij`$, $`i,j=1\mathrm{}5`$ and $`h_{ij}=h_{ji}^1`$) the corresponding 10 boundary variables (associated to thin boundary edges in Figure 1) 444Strictly speaking, the boundary connections $`h_{ij}`$ are defined as the product $`h_{ij}^{}h_{ij}^{^{\prime \prime }}`$ where $`h^{}`$ and $`h^{^{\prime \prime }}`$ are associated to half paths as follows: take the edge $`ij`$ for simplicity and assume it is oriented from $`i`$ to $`j`$. Then $`h_{ij}^{}`$ is the discrete holonomy from $`i`$ to some point in the center of the path and $`h_{ij}^{^{\prime \prime }}`$ is the holonomy from that center point to $`j`$. This splitting of variables is necessary when matching different atoms to reconstruct the simplicial amplitude. The use of this variables (wedge variables) will be crucial in Section 5. For a more detailed description of wedge variables see . and $`g_iSpin(4)`$ ($`i=1,\mathrm{},5`$) the internal connection (corresponding to the thick edges in Figure 1). According to (8) the 4-simplex BF amplitude $`4SIM_{BF}(h_{ij})`$ is given by

$$4SIM_{BF}(h_{ij})=\underset{i}{}dg_i\underset{i<j}{}\delta (g_ih_{ij}g_j).$$
(18)
With the definition of the $`B`$ fields given in (17) the constrained amplitude, $`4SIM_{const}(h_{ij})`$, formally becomes

$$4SIM_{const}(h_{ij})=\underset{i}{}dg_i\delta \left[\mathrm{Constraints}(𝒳(U_{ij}))\right]\underset{i<j}{}\delta (g_ih_{ij}g_j),$$
(19)
where $`U_{ij}=g_ih_{ij}g_j`$ is the holonomy around the triangular face (wedge) $`0ij`$ according to Figure 1. It is easy to verify, using an equation analogous to (15), that one can define the $`B`$’s by simply acting with the right invariant vector fields on the boundary connection $`h_{ij}`$. Therefore, the previous equation is equivalent to

$$4SIM_{const}(h_{ij})=\delta \left[\mathrm{Constraints}(𝒳(h_{ij}))\right]\underset{i}{}dg_i\underset{i<j}{}\delta (g_ih_{ij}g_j),$$
(20)
where we have taken the delta function out of the integral. The quantity on which the formal delta distribution acts is simply $`4SIM_{BF}(h_{ij})`$ (defined in (18)), which after integrating over the internal connection $`g_i`$, and using equation (12) becomes

$$4SIM_{BF}(h_{ij})=\underset{\rho _1\mathrm{}\rho _{10}}{}\underset{\iota _1\mathrm{}\iota _5}{}\begin{array}{c}\text{}\end{array}\begin{array}{c}\text{}\end{array},$$
(21)
where the circles represent the corresponding $`\rho `$-representation matrices evaluated on the the corresponding boundary connection $`h`$. The term on the left is a $`15j`$-symbol as in (14) while the term on the right is the trace of five intertwiners with the respective boundary connection insertions. Notice then that nodes on the two pentagonal diagrams are linked together by the value of their intertwiner.

The 4-simplex amplitude for the constraint spin foam model is then defined as the restriction of $`4SIM_{BF}(h_{ij})`$ imposed by the quantum version of the constraints (3). The latter are defined by the following set of differential equations

$$ϵ_{IJKL}𝒳^{IJ}(h_{ij})𝒳^{KL}(h_{ik})4SIM_{const}(h_{ij})=0j,k,$$
(22)
and where the index $`i=1,\mathrm{},5`$ is held fixed. The translation of the continuum constraint (3) into discrete elements associated to faces in $`\mathrm{\Delta }^{}`$ is analogous to that given in . Notice that (22) is to be thought as a condition on BF amplitudes and is not a general equation to be imposed to any 4-simplex amplitude. Recall that the strategy is to constraint the BF theory to obtain a definition of the path integral for general relativity 555We illustrate the general idea with the following simple example. Imagine that the analog of $`4SIM_{BF}`$ function (eq. (18)) is the integral
$$A=𝑑k𝑑pe^{ikx+ipy}=\delta (x)\delta (y),$$
(23)  where $`x,y[0,2\pi ]`$ represent the boundary ‘connections’. The analog of the constraint (3) is defined to be $`kp=0`$ which in turn implies the constrained amplitude to be
$$A_{const}=\delta (x+y).$$
Let us now apply the prescription used in the BF theory. We can expand the un-constrained function (23) in terms of ‘spin foam’ amplitudes
$$A=\frac{1}{4\pi ^2}\underset{n,m}{}e^{inx+imy}.$$
In this case this corresponds to Fourier expanding delta function on $`S_1\times S_1`$ (Peter-Weyl decomposition for $`U(1)\times U(1)`$). The constraint is now represented by a combination $`C`$ of right invariant vector fields on $`U(1)`$: $`C=_x_y`$. So we can now impose the constraints by means of selecting those configurations (modes) in (23) that are annihilated by $`C`$. The equation analogous to (22) is
$$(_x_y)e^{inx+imy}=(nm)e^{inx+imy}=0$$
which implies $`n=m`$ and $`A_{const}=\delta (x+y)`$. .

Equations closely related to (22) can also be obtained as the geometric restrictions on the $`B`$’s to be simple bi-vectors coming from a dual cotetrad or to characterize the geometry of a tetrahedron in 4 dimensions . In this case the equivalent of $`B`$ correspond to bivectors defined by the faces of a classical tetrahedron. Using geometric quantization one obtains the Hilbert space of states of the ‘quantum tetrahedron’ where the $`B`$’s are promoted to operators. Notice the our $`B`$ operator (15) is obtained directly from the BF path integral and one does not need to invoke any additional quantization principle. The procedure is completely analogous to the simple example of Footnote 5. A similar point of view has been taken by Reisenberger and Freidel-Krasnov in .

### 4.2 Restricted BF paths

The following procedure is very similar in spirit to the BC prescription. The essential difference is that we now require the set of restricted configurations to be contained in the set of modes of the BF amplitude, $`4SIM_{BF}`$.

There are seven equations (22) for each value of $`i=1,\mathrm{},5`$. If we consider all the equations for the 4-simplex amplitude then some of them are redundant. The total number of independent conditions is 20, in agreement with the number of classical constraints (3). For a given $`i`$ in (22) (i.e., a given tetrahedron) and for $`j=k`$ the equation becomes

$`ϵ_{IJKL}𝒳^{IJ}(h_{ij})𝒳^{KL}(h_{ij})4SIM_{const}(h_{ij})`$
$`=\left[j_{ij}^{\mathrm{}}(j_{ij}^{\mathrm{}}+1)j_{ij}^r(j_{ij}^r+1)\right]4SIM_{const}(h_{ij})=0,`$ (24)
where we have used $`\rho =j^{\mathrm{}}j^r`$ for $`j^{\mathrm{}},j^r\mathrm{Irrep}[SU(2)]`$. The previous constraints are solved by requiring the corresponding representation $`\rho _{ij}`$ to be simple, i.e., $`\rho _{ij}=j_{ij}j_{ij}`$.

This solves 10 of the 20 equations. The next non-trivial condition imposed by (22) is when $`jk`$. In this case we have

$`2ϵ_{IJKL}𝒳^{IJ}(h_{ij})𝒳^{KL}(h_{ik})4SIM_{const}(h_{ij})`$
$`=ϵ_{IJKL}\left(𝒳^{IJ}(h_{ij})+𝒳^{IJ}(h_{ik})\right)\left(𝒳^{KL}(h_{ij})+𝒳^{KL}(h_{ik})\right)4SIM_{const}(h_{ij})`$
$`=\left[\iota ^{\mathrm{}}(\iota ^{\mathrm{}}+1)\iota ^r(\iota ^r+1)\right]4SIM_{const}(h_{ij})`$
$`=0,`$ (25)
where we used the gauge invariance at the 3-valent node in the tree decomposition that pairs the representation $`\rho _{ij}`$ with the $`\rho _{ik}`$666 The gauge invariance at the node allows us to express the sum of right-invariant vector fields acting on the external ‘legs’ (see Figure 2.) as a right-invariant vector field acting on the internal representation $`\iota `$. Of course right versus left invariant vector field is a matter of convention which implies a choice of orientation. Since $`Spin(4)`$ representations are self dual we have $`𝒳_R^{IJ}(h)=𝒳_L^{IJ}(h^1)`$., and that we have already solved (4.2). In the last line we assume that the internal color of the corresponding 4-intertwiner is $`\iota =\iota ^{\mathrm{}}\iota ^r`$. This choice of tree decomposition in the case $`ij=12`$ and $`ik=13`$ is illustrated in Figure 2. The solution is clearly $`\iota =\iota \iota `$.

What happens now with any of the other two remaining conditions, for example, $`E(ij^{},ik^{})`$ for $`kk^{}`$, $`jj^{}`$ and $`j^{}k^{}`$? It seems that we ran out o possibilities of restricting the representations. Generically this equations will not be satisfied because an intertwiner that has simple $`\iota `$ in one tree decomposition has not only simple $`\iota ^{}`$’s components in a different tree decomposition and the equation would be violated. However there is a case in which this happens trivially, namely when the dimension of the invariant part of the tensor product of the four corresponding representations is unity. Let us write this condition as an equation since this corresponds to the solution of the remaining 5 independent conditions, namely

$$\mathrm{dim}\left(\mathrm{Inv}\left[\rho _{ij}\rho _{ik}\rho _{im}\rho _{ip}\right]\right)=1$$
(26)
In this case $`\iota `$ would be simple in any tree decomposition if it is simple in one particular one. Notice that this is the only solution to our constraints as a trivial consequence of the theorem proven by Reisenberger in . Our set of solutions are contained in the Barret-Crane solutions since our intertwiners agree with the BC one every time that equation (26) is satisfied. Solutions to the previous equation can be characterized as follows. Since all the $`\rho `$’s are simple we can concentrate on their right (or left) components. Assume $`j_1j_2j_3j_4`$ then the condition is $`j_1+j_2+j_3=j_4`$. Explicitly, a few examples of solutions are $`(\frac{1}{2},\frac{1}{2},2,3)`$, $`(\frac{3}{2},\frac{1}{2},\frac{5}{2},\frac{1}{2})`$, $`(1,5,226,220)`$, etc. The intertwiner color is completely determined by the face colors and a choice of tree decomposition. If we pair $`j_1`$ and $`j_2`$ in our previous example then $`\iota =j_1+j_2=j_4j_3`$. The amplitude $`4SIM_{const}`$ is independent of the tree decomposition chosen.

### 4.3 Gluing 4-simplexes

Once we have solved equations (22) for a single 4-simplex we can calculate the amplitude of any simplicial decomposition of $``$, $`\mathrm{\Delta }`$. This is achieved by putting together 4-simplexes with consistent boundary connections and gluing them together by means of integrating over the boundary data in the standard way.

Let us point out that there is a potential ambiguity in this step. We have implemented constraints in the path integral and this generally should be supplemented with the appropriate modification of the measure. This could affect the values of lower dimensional simplexes such as face and edge amplitudes. Constraints (22) act on each edge (tetrahedron) separately, heuristically one would expect a Jacobian factor to modify the edge amplitude of the model since the constraints are non linear functions of the $`B`$’s. A rigorous derivation of such factors from the path integral definition is lacking in our argument. We believe that this might shed light on the problem of the correct normalization of this type of spin foam models. A more detailed study of this issue is being explored in .

If we do so then we end up with

$`Z_{const}(\mathrm{\Delta })={\displaystyle \underset{𝒞_f:\{f\}\rho _f^s}{}}{\displaystyle \underset{f\mathrm{\Delta }^{}}{}}\mathrm{\Delta }_\rho {\displaystyle \underset{e\mathrm{\Delta }^{}}{}}A_e{\displaystyle \underset{v\mathrm{\Delta }^{}}{}}\begin{array}{c}\text{}\end{array},`$ (28)
where $`\rho _f^s`$ denote the set of representations selected by conditions (22), $`\iota _e^s`$ are the corresponding colors of intertwiners, and $`A_e`$ is the appropriate edge amplitude (undetermined in our prescription).

### 4.4 Volume

In this section we discuss a rather puzzling feature of the model we have defined above.

If we consider boundaries, then the spin-network states induced as boundaries of spin foams are four-valent and the representations of the corresponding edges satisfy (26). Now using the standard definition of the volume operator on this set of states we obtain an identically zero result, i.e, the 3-volume operator $`V_{(3)}`$ annihilates the states that solve (22). The reason is that the volume is given by

$$V_{(3)}^2[(𝒳_i+𝒳_j)^2,(𝒳_i+𝒳_k)^2],$$
(29)
where the square is taken using the internal metric $`\delta _{IJ}`$. The solutions to the constraints happen to diagonalize both operators in the commutator which implies $`V_{(3)}=0`$.

## 5 Degenerate sector

As shown in , constraints (3) correspond to the non-degenerate phase of solutions of the general constraints (i.e., phase with $`e0`$). In  Reisenberger explicitly solved the constraints in the degenerate sectors and showed that, in these cases, the action reduces to

$$S_{deg}^\pm =B_i^r(F_i(A^r)\pm V_i^jF_j(A^{\mathrm{}})),$$
(30)
where the upper index $`r`$ (respectively $`\mathrm{}`$) denotes the self-dual (respectively anti-self-dual) part of $`B`$ and $`A`$ in the internal space, and $`VSO(3)`$.

Let us concentrate in the sector with the minus sign in the previous expression. Then it is straightforward to define the discretized path integral along the same lines as BF theory in Section 3. The result is

$$𝒵(\mathrm{\Delta })=\underset{f\mathrm{\Delta }^{}}{}dB_f^{r\left(3\right)}dv_f\underset{e\mathrm{\Delta }^{}}{}dg_e^{\mathrm{}}dg_e^re^{i\mathrm{Tr}\left[B_f^rU_f^rv_fU_f^\mathrm{}1v_f^1\right]}.$$
(31)
Integrating over the $`B`$ field we obtain

$$Z(\mathrm{\Delta })=\underset{e\mathrm{\Delta }^{}}{}dg_e^{\mathrm{}}dg_e^r\underset{f\mathrm{\Delta }^{}}{}dv_f\delta ^{(3)}(g_{e_1}^r\mathrm{}g_{e_n}^rv_f(g_{e_1}^{\mathrm{}}\mathrm{}g_{e_n}^{\mathrm{}})^1v_f^1),$$
(32)
where $`dg_e^{\mathrm{}}`$, $`dg_e^r`$, and $`dv_f`$ are defined in terms of the $`SU(2)`$ Haar measure and the delta function $`\delta ^{(3)}`$ denotes an $`SU(2)`$ distribution.

In order to obtain the corresponding state-sum it is easier to concentrate on a single 4-simplex amplitude. Furthermore, we start by the wedge shown in Figure 3. In this figure we represent one of the 10 wedges that form a 4-simplex atom (see Figure 1). Both the internal connection $`g_{ij}`$ ($`g_{ij}=g_{ji}^1`$) and the boundary connection variables $`h_{ij}`$ ($`h_{ij}=h_{ji}^1`$) are in $`Spin(4)`$, while $`u_{ljki}SU(2)Spin(4)`$ is an auxiliary variable. The $`SU(2)`$ subgroup is defined as the diagonal insertion $`ug=(ug^{\mathrm{}},ug^r)`$. The wedge amplitude is defined as

$$w=𝑑u_{ljki}\delta ^{(6)}(g_{ki}h_{il}u_{ljki}h_{lj}g_{jk})$$
(33)
according to the notation in Figure 3 and where the $`\delta ^{(6)}`$ denotes a $`Spin(4)`$ delta distribution.

Any face in the 2-complex will be defined by as many such wedges as 4-simplexes share the corresponding face. Figure 4 illustrates the case for a triangular face. The vertices $`1,2`$ and $`3`$ correspond to the centers of the three 4-simplexes sharing the face. The dotted line denotes the region along which the boundary of the three atoms (Figure 1) join.

It is easy to check that integrating over all but one boundary variables $`h_{ij}`$, the contribution of a combination of wedges forming a face $`f\mathrm{\Delta }^{}`$ is given by

$$𝑑u_f𝑑h\delta ^{(6)}(U_fhu_fh^1),$$
(34)
where $`U_fSpin(4)`$ is the discrete holonomy, $`u_fSU(2)Spin(4)`$ is a product of the $`u_w`$ associated to the corresponding wedges and $`h`$ is the remaining boundary connection. In the case shown in Figure 4, $`U_f=g_{1i}g_{i2}g_{2j}g_{j3}g_{3k}g_{k1}`$, $`u_f=u_3u_2u_1`$, and $`h=h_{kl}`$. Using that $`Spin(4)=SU(2)\times SU(2)`$ and the definition of the $`SU(2)`$ subgroup where $`u`$ lives, the integral over $`u_f`$ of the previous equation becomes

$`{\displaystyle 𝑑u_f𝑑h^{\mathrm{}}𝑑h^r\delta ^{(3)}(U_f^{\mathrm{}}h^{\mathrm{}}u_fh^\mathrm{}1)\delta ^{(3)}(U_f^rh^ru_fh^{r1})}`$
$`={\displaystyle d(h^{\mathrm{}}h^{r1})\delta ^{(3)}(U_f^{\mathrm{}}h^{\mathrm{}}h^{r1}U_f^{r1}h^rh^\mathrm{}1)},`$ (35)
where we have used that $`dh=dh^{\mathrm{}}dh^r`$ and $`\delta ^{(6)}(g)=\delta ^{(3)}(g^{\mathrm{}})\delta ^{(3)}(g^r)`$ as well as the invariance and normalization of the Haar measure. The previous face amplitude coincides with that in (32) if we define $`v_f=h^{\mathrm{}}h^{r1}`$. Therefore, (33) defines the wedge amplitude of (32).

Now we can write the analog of equation (21) for the 4-simplex amplitude, $`4SIM_{Deg}(h_{ij})`$, putting together the 10 corresponding wedges and integrating over the internal $`g`$’s, namely

$$4SIM_{Deg}(h_{ij})=\underset{\rho _1\mathrm{}\rho _{10}}{}\underset{\iota _1\mathrm{}\iota _5}{}\begin{array}{c}\text{}\end{array}\begin{array}{c}\text{}\end{array},$$
(36)
where the dark dots denote integration over the $`SU(2)`$ diagonal subgroup, and the white circles represent the boundary connections. To keep the diagrammatic notation simple we have dropped some labels. The next step is to perform the integration over the $`u`$’s. We concentrate on a single intertwiner in (36), i.e., a single node in the pentagonal diagram on the right of the previous equation. Using the orthogonality of $`SU(2)`$ unitary irreducible representations and the fact that the representations $`\rho `$ of $`Spin(4)`$ are of the form $`\rho =jk`$ for $`j,k`$ $`SU(2)`$ unitary irreducible representations we have

$$\begin{array}{c}\text{}\end{array}=\frac{\delta _{j_1,k_1}\mathrm{}\delta _{j_4,k_4}}{(2j_1+1)\mathrm{}(2j_4+1)}\begin{array}{c}\text{}\end{array},$$
(37)
where we represent representations $`\rho =jk`$ as parallel lines, the symbol $`R`$, $`L`$ in the circles on the left denotes $`h^r`$, $`h^{\mathrm{}}`$ $`SU(2)`$-representation matrices, $`RL`$ on the right denotes the product $`h^rh^1\mathrm{}`$, and the dark dots are subgroup integrations.

The Kronecker deltas in the previous equation implies that the $`\rho `$’s labeling faces must be simple, i.e., $`\rho =jj`$.

Finally, it is easy to verify that when gluing various 4-simplex atoms together by means of integrating over matching boundary connections the integration simply set $`\iota _R=\iota _R^{}`$, $`\iota _L=\iota _L^{}`$ where $`\iota `$, $`\iota ^{}`$ are the intertwiners corresponding to the tetrahedron shared by the two 4-simplexes. Notice that no simplicity condition is imposed on $`\iota `$. Now consider an arbitrary face bounded by $`n`$ edges. Such a face is made up $`n`$ wedges. Therefore there is a factor $`\mathrm{\Delta }_{jj}^{2n}=(2j+1)^{2n}`$ coming from the delta function mode expansion in (recall Footnote (2).), a factor $`(2j+1)^n`$ from the factors in (37), and finally a factor $`(2j+1)^n`$ from the boundary connection integrations in the gluing. This results in a face amplitude equal to unity.

Putting all this together one gets a spin foam model were only face representations are constrained to be simple while intertwiners are arbitrary. Explicitly

$`Z_{deg}(\mathrm{\Delta })={\displaystyle \underset{𝒞_f:\{f\}\rho _f^s}{}}{\displaystyle \underset{𝒞_e:\{e\}\{\iota _e\}}{}}{\displaystyle \underset{v\mathrm{\Delta }^{}}{}}\begin{array}{c}\text{}\end{array}.`$ (39)
This is precisely the spin foam obtained in ! This model was obtained as a natural modification of the GFT that defines a variant of the BC model. Here we have rediscovered the model from the systematic quantization of $`S_{deg}^{}`$ defined in (30). This establishes the relation of the model with a classical action! It corresponds to spin foam quantization of the ‘$``$’ degenerate sector of $`SO(4)`$ Plebanski’s theory.

The $`+`$ sector action (30) can be treated in a similar way. The only modification is that of the subgroup. Instead of using the diagonal insertion defined above one has to define $`uSU(2)Spin(4)`$ so that $`ug=(ug^{\mathrm{}},u^1g^r)`$.

We have restricted to simplicial decompositions but all this should be generalizable along the lines of reference  for arbitrary cellular decompositions of $``$. This generalization seems straightforward although it should be investigated in detail.

To conclude this section let us notice that the allowed 4-simplex configurations of the model of Section 4 are fully contained in the set of 4-simplex configurations of the model obtained here. We come back to this issue in the following section.

## 6 Discussion

The principal idea behind this work was to study the spin foam quantization of Plebanski formulation of gravity by restricting the paths that appear in the $`SO(4)`$ BF theory. This strategy is supported by the fact that Plebanski’s action can be thought of as the $`SO(4)`$ BF theory, supplemented by certain constraints on the B field. Gravity in the Palatini formulation is obtained as one of the non-degenerate sectors of the solutions to the classical constraints. In the model introduced in section 4 we defined a prescription for implementing these constraints by restricting the set of histories of the BF theory to those satisfying the ‘quantum analog’ of (3). Solution configurations of a single 4-simplex in the model are special 4-simplex BF configurations. Even though the 4-simplex configurations appearing here are a sub-set of the Barrett-Crane configurations, the great majority of the Barrett-Crane configurations (independently of the normalization chosen) are excluded by the requirement that they be BF configurations. The nature of the constraints in the BC model is essentially the same as the ones defined here. The difference is that in the BC case constraints are implemented on a single intertwiner while here we keep track of the fact that intertwiners appear in pairs in the BF state sum (see (12)).

There are alternative ways to motivate the definition of the Barrett-Crane model which are independent of the line of thought used here. There is also evidence that relate it to a theory of quantum gravity. However, we believe that this work shows that there is no obvious way to interpret it as the quantization of Plebanski’s action. Using reasonable definitions we have shown that one obtains a more restrictive state sum.

In the context of the BC model, reference  shows how one can restrict the states of the ‘quantum tetrahedron’ so that fake tetrahedra are ruled out of the state sum. In our context this amounts to resolving the ambiguity between the $`ee`$ and $`{}_{}{}^{}(ee)`$ solutions of the constraints (see (4)) at the quantum level. In  it is shown that the ‘correct’ configurations are selected by imposing the so-called chirality constraint which is automatically satisfied at the quantum level because it can be written as the commutator of the simplicity constraints (22). It is also shown that the two spaces of solutions ($`ee`$ and $`{}_{}{}^{}(ee)`$ respectively) intersect on the set of configurations for which $`V_{(3)}=0`$. It is easy to see that all this can be translated to our context. The vanishing of the volume operator implies that in our model one can not distinguish the two type of configurations and that the ambiguity (4) remains at the quantum level.

The model of Section 4 contains only degenerate configurations in the sense that spin-network states on the boundaries have zero volume. This shows that the model cannot reproduce any of the semi-classical states of general relativity. Somehow our definition of the constraints at the level of the state-sum are so strong that non-degenerate configurations have been eliminated. Some evidence supporting this view can be obtained considering the following argument. Constraints are implemented locally on each 4-simplex; therefore, we can concentrate on a single 4-simplex to analyze their action. If we do so, then we conclude that all the 4-simplex configurations in the model (28) are entirely contained in the set of 4-simplex configurations of the model found in Section 5.

Is there a way out? If we maintain the point of view of defining the model starting from Plebanski’s action then the problem can be traced back to our definition of constraints. As it was pointed out in  there are two ways to write Plebanski’s constraints in the non-degenerate sector (i.e., when $`e0`$ in (3)). If we stick with (3) then one can try to change the definition of the $`B`$ operators. One possibility would be to change right-invariant by left invariant vector fields in the definition (17). One can do this consistently only if one changes the orientation of the 4-simplex in which case the final result remains the same. One can try to use the sum of right-invariant and left invariant vector fields. This is certainly a possibility, and actually converges faster to the value of $`B`$ in the sense of equation (15) when $`U1`$. However, since the constraints are quadratic in the $`B`$ there will be cross-terms in (22). This terms cannot be expressed in terms of $`Spin(4)`$ Casimir operators, and consequently, the constraints cannot be solved in terms of simple restrictions on the set of representations involved in the state-sum. There seem to be no obvious way to use (3) and avoid the discouraging results of Section 4.

The other possibility is to use

$$ϵ^{\mu \nu \rho \sigma }B_{\mu \nu }^{IJ}B_{\rho \sigma }^{KL}ϵ^{IJKL},$$
(40)
where the difference with (3) is in the fact that we have traded internal with space-time indices. However with this choice, the connection with the BC model becomes much more uncertain. Notice that in this form, the constraints have free Lie-algebra indices and therefore cannot be written as Casimir operators as they stand. This version of the constraints has been studied in the literature. Such constraints have been incorporated in a spin foam model of Riemannian general relativity in terms of self-dual variables by Reisenberger in . In the context of $`SO(4)`$ Plebanski’s action a model along this lines has been defined by Freidel et al. in . But all these models are quite different from the BC model. We believe that this shows that there is no obvious means of interpreting the BC model as a spin foam quantization of Plebanski’s theory.

Let us conclude by analyzing the results of the last section. In Section 5 we quantized the degenerate sectors of Plebanski’s action in a fairly straightforward way. In this case we do not impose any constraints and the state sum follows directly from the discretized definition of the path integral of the theory. There are no ambiguities in lower dimensional simplex-amplitudes. The model turns out to be precisely the one introduced by De Pietri, Freidel, Krasnov and Rovelli in . This work establishes a clear connection between that model and the effective action corresponding to one of the degenerate sectors of Plebanski’s action.

Finally the model is well defined, is not topological and has a clear connection to a continuous action. It is somehow between the theory we want to define and the simpler theories we understand well but do not have local excitations (such as BF theory and gravity in lower dimensions). From this viewpoint we believe that it might be useful to explore its properties as a ‘toy model’ for understanding open issues in the spin foam approach to quantum gravity. Among these is the very important problem of the continuum limit (i.e., the issue of summing-over versus refining discretizations) and the interpretation of the path integral in the diffeomorphism invariant context (time evolution versus the projector/extractor operator on physical states).

## 7 Acknowledgments

I thank Abhay Ashtekar, Martin Bojowald and Amit Ghosh for stimulating discussion and insightful questions. I also thank John Baez, Rodolfo Gambini, Jorge Pullin, Michael Reisenberger and Carlo Rovelli for motivating remarks. I am grateful to Abhay for many suggestions that have improve the presentation of this manuscript. This work was supported in part by NSF Grant PHY-0090091, and Eberly Research Funds of Penn State.


text:
# Transition to Fulde-Ferrel-Larkin-Ovchinnikov phases near the tricritical point : an analytical study

## I INTRODUCTION

Although they have been proposed a long time ago, Fulde-Ferrel- Larkin-Ovchinnikov (FFLO) phases  are still the subject of a continuing interest. Indeed the existence of these phases is a fairly remarkable phenomenon since they correspond to a spontaneous symmetry breaking of the standard BCS superfluid phase in the presence of an effective field, inducing a difference in chemical potential between the two populations involved in the formation of Cooper pairs. This symmetry breaking leads to an inhomogeneous superfluid with a space dependent order parameter, while the applied field is perfectly homogeneous. This situation is analogous to the appearance of vorticity in type II superconductors, but in this latter case the effect is due to the coupling of the field to particle currents while in FFLO phases only the coupling to the spin of the pairing fermions is involved. In standard superconductors the coupling to the orbital degrees of freedom is much stronger than the coupling to the spins. Hence the upper critical field is due to the orbital coupling and the FFLO phases can not be observed, since they should appear at much higher field. However in heavy fermions superconductors the strength of these two couplings is comparable, which could make possible the observation of FFLO phases. Nevertheless their sensitivity to impurities could be a major problem. Another possible direction to eliminate the orbital coupling is to consider lower dimensional superconductors, in a geometry where the currents would have to flow in an actually prohibited direction. Organic compounds or cuprate superconductors are interesting systems in this respect. And indeed very recently the FFLO state has been claimed to be observed in a quasi-two-dimensional organic compound  . On the other hand earlier possible observations in heavy fermion compounds  have not been undisputed. We note in particular that the analysis of experimental results relies very often heavily on the theoretical results, but we will see that the situation is not completely satisfactory in this respect.

Another class of physical systems where FFLO phases could be observed is coming up quite recently. These are the ultracold fermionic gases. As it is well known remarkably low temperatures have been obtained on bosonic gases, leading in particular to the observation of Bose-Einstein condensation in alkali ultracold gases. More recently fermionic gases have been cooled down in the degenerate regime  and reaching a BCS superfluid transition in these systems seems a reasonable possibility  . However in the systems considered for observing this transition there is, in contrast to electronic spin relaxation in superconductors, no fast relaxation mechanism to equalize the populations of the two fermion species involved in the formation of Cooper pairs. Hence one should have no limitation to the effective field in these systems since the number of atoms in the populations can be in principle obtained at will. So the difference in atomic populations looks as a very promising control parameter. On the other hand if this parameter is not fully controlled this might very well be a major difficulty in reaching the BCS transition in these systems  .

In contrast with this raising experimental interest there are still theoretical problems with the precise nature of the possible phases. Specifically the basic FFLO instability corrresponds to have pairs formed with a total nonzero momentum $`𝐪_0`$ instead of forming pairs ($`𝐤,𝐤`$) with zero total momentum as in the standard BCS phase. This gives rise to a spatial dependence $`\mathrm{exp}(i𝐪_0.𝐫)`$ for the order parameter, which leaves a degeneracy with respect to the orientation of $`𝐪_0`$. This has been investigated by  Larkin and Ovchinnikov (LO) who looked how it is lifted right below the critical field. In this case, when considering the spatial dependence of the order parameter $`\mathrm{\Delta }(𝐫)`$, one can restrict the investigation to the subspace generated by linear combination of the plane waves $`\mathrm{exp}(i𝐪_0.𝐫)`$ with all possible directions $`𝐪_0`$. At T = 0, LO looked for periodic structures and found that the energetically favored result is a second order transition to a one-dimensional ’planar’ texture $`\mathrm{\Delta }(𝐫)\mathrm{cos}(𝐪_0.𝐫)`$. However they left open in their paper the possibility of a first order transition. Actually when considering the three-dimensional ’cubic’ texture $`\mathrm{\Delta }(𝐫)\mathrm{cos}(q_0x)+\mathrm{cos}(q_0y)+\mathrm{cos}(q_0z)`$ they found that it is energetically unfavorable compared to the normal state. Nevertheless they obtained from the gap equation that a nonzero solution for this order parameter exists *above* the FFLO transition line. In terms of the expansion of the free energy in powers of the order parameter, schematically $`\mathrm{\Omega }=\alpha _2\mathrm{\Delta }^2+\alpha _4\mathrm{\Delta }^4`$, this situation corresponds to a positive coefficient $`\alpha _2`$ for the second order term and a negative coefficient $`\alpha _4`$ for the fourth order one, just the opposite of the standard Landau-Ginzburg expansion below the transition. The LO evaluation of the related free energy corresponds actually to the maximum $`\alpha _2^2/4|\alpha _4|`$ of this free energy. Beyond this maximum the free energy decreases, and it would go to $`\mathrm{}`$ if one would consider only the second and fourth order terms. Naturally one has to include the effect of all higher order terms in order to find the value of the free energy for large values of the order parameter. However, at the FFLO transition line and slightly above it, the free energy becomes negative for values of the order parameter (specifically for $`\mathrm{\Delta }^2=\alpha _2/|\alpha _4|`$) where only the second and fourth order terms have to be kept in the expansion. This shows definitely that the transition toward the FFLO phase occurs *above* the standard second order FFLO transition line and is actually first order, because one can display in this range a solution which has a lower free energy than the normal state. On the other hand in order to obtain consistently the order parameter which gives the lowest free energy one has naturally to take into account higher order terms. In particular it is by no means obvious that the cubic phase is the stable one. Moreover when higher order terms are considered there are no reason anymore to restrict the search to the subspace generated by the plane waves $`\mathrm{exp}(i𝐪_0.𝐫)`$.

In order to explore more fully this difficult problem with easier conditions, it is better to be able to proceed to some kind of expansion. This can be done if, instead of working at $`T=0`$, one explores the vicinity of the tricritical point (TCP), where the FFLO transition line starts. It is located at $`T_{\mathrm{tcp}}/T_{c0}=0.561`$ where $`T_{c0}`$ is the critical temperature for $`\overline{\mu }=0`$ , with $`2\overline{\mu }=\mu _{}\mu _{}`$ being half the chemical potential difference between the two fermionic populations forming pairs. The corresponding effective field is $`\overline{\mu }_{\mathrm{tcp}}/T_{c0}=1.073`$. At this point, in the free energy expansion, both the coefficient of the second ordre term $`\alpha _2`$ and the coefficient of the fourth ordre term $`\alpha _4`$ vanish. Naturally $`\alpha _2`$ is zero just because the TCP is on the standard second order phase transition line. On the other hand following this transition line one has $`\alpha _4>0`$ for $`T>T_{\mathrm{tcp}}`$ and $`\alpha _4<0`$ for $`T<T_{\mathrm{tcp}}`$. The change of sign of $`\alpha _4`$ at $`T_{\mathrm{tcp}}`$ is the origin of the FFLO instability. Clearly, by continuity, the equilibrium order parameter will be small in the vicinity of the TCP since we know that it is zero just above the TCP on the second order line and it is small right below it in the superfluid phase. Therefore a power expansion of the free energy will be enough to find it. On the other hand, since the second and fourth order terms are zero at the TCP, we have clearly to expand at least up to sixth order, but this will prove to be enough because the corresponding coefficient is positive and not small. Similarly the optimum wavevector $`𝐪_0`$ corresponding to the FFLO phase will be small in the vicinity of the TCP since it is zero just above it on the second order line. This allows to proceed to a gradient expansion of the free energy. Since only even powers of the wavevector can enter and we look for a minimum as a function of this wavevector, we have to expand at least up to fourth order in gradient, but this will prove again to be enough.

This line of thought has actually already been followed by Houzet et al. , who have performed this expansion for the free energy and explored the result numerically. They have found that the energetically favored phase is the one-dimensional planar order parameter found by LO at $`T=0`$, but that the transition is actually slightly first order, instead of second order as found by LO at $`T=0`$. Our purpose in the present paper is rather to proceed to an analytical study of this problem. Indeed there are infinitely many possible order parameters in competition. And our aim, in considering the vicinity of the TCP, is to find the important ingredients which are responsible for the selection of the actual stable state and obtain a better physical understanding, having in particular in mind the generalization to more complicated situations. Hence our paper is complementary to their work. In particular we obtain a first order transition to the one- dimensional planar order parameter, but we will be able to analyze the reasons which favor this phase. The transition to this planar order parameter has been actually explored numerically down to $`T=0`$ by Matsuo et al.  . They have used quasiclassical equations and found that the transition keeps first order down to low temperature, but eventually goes to second order in agreement with LO. On the other hand since we know that at $`T=0`$ the cubic phase is more stable than the planar one, the question of the stablest phase at low temperature is still unsolved. Finally this first order transition to the planar phase in the three-dimensional case is in contrast with the results of Burkhardt and Rainer  who found it to be second order in a two-dimensional space.

In the following section, for completeness and to set up our notations, we rederive the expression  of the free energy. After considering in section III some simple situations, we explore in details the minimization of this free energy. This is done in section IV by restricting our search to the LO subspace for the order parameter, which is made of arbitrary superpositions of plane waves. In section V we show that our results are only slightly modified when we release this restriction. Throughout the paper we restrict ourselves to the simplest BCS scheme, namely we will consider the free energy corresponding to a weak coupling isotropic Fermi system, ignoring in particular any Fermi liquid effect. Moreover we concentrate on the three-dimensional case which leads to a first order transition, and make only occasionally comparison with the two-dimensional situation where the transition is second order.

## II THE FREE ENERGY

There are various ways to obtain the explicit expression for the free energy we need in the vicinity of the TCP . In practice it is convenient to use the fact  that, by varying the free energy with respect to $`\mathrm{\Delta }(𝐫)`$ one finds the gap equation, which is also easily obtained from Gorkov s equations, as done for example by LO  . The integral form of these equations is  , with standard notations:

$`G_{}(𝐫,𝐫^{})=G_{}^0(𝐫𝐫^{}){\displaystyle 𝑑𝐫_1G_{}^0(𝐫𝐫_1)\mathrm{\Delta }(𝐫_1)F^+(𝐫_1,𝐫^{})}`$ (1)
$`F^+(𝐫,𝐫^{})={\displaystyle 𝑑𝐫_1\overline{G_{}^0}(𝐫𝐫_1)\mathrm{\Delta }^{}(𝐫_1)G_{}(𝐫_1,𝐫^{})}`$ (2)
with, for the Fourier transforms of the free fermions thermal propagators, $`G_{}^0(𝐤)=(i\omega _n\xi _𝐤+\overline{\mu })^1`$ and $`\overline{G_{}^0}(𝐤)=(i\omega _n\xi _k\overline{\mu })^1`$ where $`\xi _𝐤`$ is the kinetic energy measured from the Fermi surface for $`\overline{\mu }=0`$ and $`\omega _n=\pi T(2n+1)`$ are Matsubara frequencies. The order parameter is given by the self- consistency relation :

$`\mathrm{\Delta }^{}(𝐫)=VT{\displaystyle \underset{n}{}}F^+(𝐫,𝐫)`$ (3)
We expand Eq. (1-2) up to fifth order in $`\mathrm{\Delta }(𝐫)`$. We introduce the Fourier transform $`\mathrm{\Delta }_𝐪=d𝐫\mathrm{\Delta }(𝐫)\mathrm{exp}(i𝐪.𝐫)`$ of the order parameter. As explained in the introduction we proceed also to an expansion in the wavevector $`𝐪`$ of the order parameter since we know that its relevant values will be small in the vicinity of the TCP. More precisely we will see that, in order to obtain a coherent expansion, it is enough to go only up to fifth order terms in overall power of $`\mathrm{\Delta }`$ and $`𝐪`$. This means that, in the gap equation, we have to expand the first order term in $`\mathrm{\Delta }`$ only to fourth order in $`𝐪`$. Similarly the third order term in $`\mathrm{\Delta }`$ has to be expanded only to second order in $`𝐪`$ and the fifth order term in $`\mathrm{\Delta }`$ can be calculated to zeroth order in $`𝐪`$. For example in order to find the third order term, we have to expand up to second order in wavevectors :

$`{\displaystyle \underset{k}{}}\overline{G^0}(𝐤)G^0(𝐤+𝐪_1)\overline{G^0}(𝐤+𝐪𝐪_3)G^0(𝐤+𝐪)\mathrm{\Delta }_{𝐪_1}^{}\mathrm{\Delta }_{𝐪_2}\mathrm{\Delta }_{𝐪_3}^{}`$ (4)
where we have used $`𝐪_1+𝐪_3=𝐪+𝐪_2`$ and we have omitted the unnecessary spin index. In the expansion appear the following numerical coefficients:

$`a_0(\overline{\mu },T)={\displaystyle \frac{1}{N_0V}}2\pi T\mathrm{}[{\displaystyle \underset{n=0}{}}{\displaystyle \frac{1}{\overline{\omega }_n}}]`$ (5)
$`a_2(\overline{\mu }/T)=\overline{\mu }^22\pi T\mathrm{}[{\displaystyle \underset{n=0}{}}{\displaystyle \frac{1}{\overline{\omega }_n^3}}]`$ (6)
$`a_4(\overline{\mu }/T)=\overline{\mu }^42\pi T\mathrm{}[{\displaystyle \underset{n=0}{}}{\displaystyle \frac{1}{\overline{\omega }_n^5}}]`$ (7)
where $`\overline{\omega }_n=\omega _ni\overline{\mu }`$, and the summation for $`a_0`$ has to be cut-off in the standard BCS way. The simple second order transition line to a standard BCS superfluid with space independent order parameter is given by $`a_0(\overline{\mu },T)=0`$. Below the TCP it corresponds to a spinodal transition line, at which the normal state becomes absolutely unstable against a transition toward a space independent order parameter. The domain $`a_0(\overline{\mu },T)>0`$ corresponds to the region of the $`(\overline{\mu },T)`$ phase diagram above this line, and it is the domain where we will look for other transitions. In practice we can see $`a_0(\overline{\mu },T)`$ as a measure of the distance from the spinodal line in the $`(\overline{\mu },T)`$ plane. Explicitely if we define $`T_{sp}(\overline{\mu }/T)`$ the spinodal temperature as a function of the ratio $`\overline{\mu }/T`$, we have $`a_0(\overline{\mu },T)=\mathrm{ln}[T/T_{sp}(\overline{\mu }/T)]`$. We will not need to explicit further this distance. As indicated in the introduction we have by definition $`a_2(\overline{\mu }/T)=0`$ at the TCP and it is small in the vicinity of this point. For $`(\overline{\mu }/T)>(\overline{\mu }/T)_{\mathrm{tcp}}=1.91`$, we have $`a_2>0`$ and $`a_2<0`$ for $`(\overline{\mu }/T)<(\overline{\mu }/T)_{\mathrm{tcp}}`$. Finally $`a_4(\overline{\mu }/T)=0.114`$ at the TCP (while it is negative near $`\overline{\mu }=0`$ and goes to $`0.25`$ when $`T0`$). With these notations the gap equation in the vicinity of the TCP reads :

$`\mathrm{\Delta }_𝐪[a_0{\displaystyle \frac{1}{3}}a_2Q^2+{\displaystyle \frac{1}{5}}a_4Q^4]{\displaystyle \underset{𝐪_i}{}}\mathrm{\Delta }_{𝐪_1}\mathrm{\Delta }_{𝐪_2}^{}\mathrm{\Delta }_{𝐪_3}[{\displaystyle \frac{1}{2}}a_2{\displaystyle \frac{1}{6}}a_4(Q^2+2Q_1^2+2Q_3^2𝐐.(𝐐_1+𝐐_3)+3𝐐_1.𝐐_3)]`$ (8)
$`+{\displaystyle \frac{3}{8}}a_4{\displaystyle \underset{𝐪_i}{}}\mathrm{\Delta }_{𝐪_1}\mathrm{\Delta }_{𝐪_2}^{}\mathrm{\Delta }_{𝐪_3}\mathrm{\Delta }_{𝐪_4}^{}\mathrm{\Delta }_{𝐪_5}=0`$ (9)
where we have used the dimensionless wavevector $`𝐐=𝐪v_F/2\overline{\mu }`$ and expressed $`\mathrm{\Delta }_𝐪`$ in units of $`\overline{\mu }`$. Also the momentum conservation is assumed in the summations, that is $`𝐪_1+𝐪_3=𝐪+𝐪_2`$ in the third order term and $`𝐪_1+𝐪_3+𝐪_5=𝐪+𝐪_2+𝐪_4`$ in the fifth order one. The above expression can be checked against the case of the simple Fulde-Ferrell state $`\mathrm{\Delta }(𝐫)=\mathrm{exp}(i𝐪_0.𝐫)`$ where a single wavevector enters.

Now the above gap equation is obtained by minimizing the following free energy difference $`\mathrm{\Omega }`$ between the superfluid and the normal state :

$`\mathrm{\Omega }={\displaystyle \underset{𝐪}{}}|\mathrm{\Delta }_𝐪|^2[a_0{\displaystyle \frac{1}{3}}a_2Q^2+{\displaystyle \frac{1}{5}}a_4Q^4]{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐪_i}{}}\mathrm{\Delta }_{𝐪_1}\mathrm{\Delta }_{𝐪_2}^{}\mathrm{\Delta }_{𝐪_3}\mathrm{\Delta }_{𝐪_4}^{}[{\displaystyle \frac{1}{2}}a_2{\displaystyle \frac{1}{6}}a_4(Q_1^2+Q_2^2+Q_3^2+Q_4^2+𝐐_1.𝐐_3+𝐐_2.𝐐_4)]`$ (10)
$`+{\displaystyle \frac{1}{8}}a_4{\displaystyle \underset{𝐪_i}{}}\mathrm{\Delta }_{𝐪_1}\mathrm{\Delta }_{𝐪_2}^{}\mathrm{\Delta }_{𝐪_3}\mathrm{\Delta }_{𝐪_4}^{}\mathrm{\Delta }_{𝐪_5}\mathrm{\Delta }_{𝐪_6}^{}`$ (11)
where we have the momentum conservation $`𝐪_1+𝐪_3=𝐪_2+𝐪_4`$ in the fourth order term while $`𝐪_1+𝐪_3+𝐪_5=𝐪_2+𝐪_4+𝐪_6`$ holds in the sixth order one. We have used symmetry and momentum conservation to present the fourth order term in a symmetrical way. This expression Eq.(7) is just the free energy we were looking for. We have considered here the 3-D case. For a two- dimensional system the angular averages found in the calculation are different. The result is simply obtained from the above one by multiplying the $`Q^2`$ terms by $`3/2`$ and the $`Q^4`$ terms by $`15/8`$.

## III SIMPLE CASES

Let us first consider some simple situations. If we consider an homogeneous order parameter, that is $`𝐪=0`$, we have merely $`\mathrm{\Omega }=a_0\mathrm{\Delta }^2a_2\mathrm{\Delta }^4/4+a_4\mathrm{\Delta }^6/8`$ . If we want to have this free energy negative for $`a_0>0`$ we need to have $`a_2>0`$, that is to be at temperature below the TCP. In this case $`\mathrm{\Omega }>0`$ when $`a_0>a_2^2/8a_4`$ and we reach a first order transition for $`a_0=a_2^2/8a_4`$, with a non zero order parameter $`\mathrm{\Delta }^2=a_2/a_4`$. This is the standard first order Pauli limiting transition. We consider next the possibility of a second order transition. In this case only the second order term in Eq.(7) is relevant. The location of this transition is given by $`a_0=\frac{1}{3}a_2Q^2\frac{1}{5}a_4Q^4`$. Below the TCP, where $`a_2>0`$, we can find $`a_0>0`$ for non zero wavevector $`𝐐`$, that is we will find an FFLO phase. Precisely the optimal wavevector is $`Q_0^2=\frac{5}{6}a_2/a_4`$ and the corresponding maximal $`a_0`$ is $`a_0=\frac{5}{36}a_2^2/a_4`$. We see that this value is larger than the one we just found for the standard Pauli limiting transition. Thus as expected the FFLO transition happens first and overtakes the first order transition. Finally it is natural and interesting to try to generalize the two above situations and consider the possibility of a first order transition for an order parameter with a single wavector component $`\mathrm{\Delta }_𝐪`$. With the shorthand $`\mathrm{\Delta }_𝐪\mathrm{\Delta }`$ the free energy writes :

$`\mathrm{\Omega }=[a_0{\displaystyle \frac{1}{3}}a_2Q^2+{\displaystyle \frac{1}{5}}a_4Q^4]\mathrm{\Delta }^2{\displaystyle \frac{1}{4}}[a_22a_4Q^2]\mathrm{\Delta }^4+{\displaystyle \frac{1}{8}}a_4\mathrm{\Delta }^6`$ (12)
Minimizing first with respect to $`Q^2`$ we obtain for the extremum the condition $`Q^2=\frac{5}{6}a_2/a_4\frac{5}{4}\mathrm{\Delta }^2`$, which implies that must have $`\mathrm{\Delta }^2\frac{2}{3}a_2/a_4`$ otherwise we are back to the homogeneous situation and the Pauli limiting transition. Inserting this value for $`Q^2`$ in Eq.(8) we find $`\mathrm{\Omega }=[a_0\frac{5}{36}a_2^2/a_4]\mathrm{\Delta }^2+\frac{1}{6}a_2\mathrm{\Delta }^4\frac{3}{16}a_4\mathrm{\Delta }^6`$. We are naturally interested in finding a transition higher than the standard FFLO. This means we are looking for $`a_0>\frac{5}{36}a_2^2/a_4`$, so the first term in $`\mathrm{\Omega }`$ is positive. But one sees that the sum of the last two terms is also positive for $`\mathrm{\Delta }^2\frac{2}{3}a_2/a_4`$. Therefore we have not been able to improve the standard FFLO solution. However we have clearly not done our best in this direction.

Before trying to improve in this way, it is convenient to simplify our expression for the free energy by taking reduced units for the order parameter and the wavector, which come out naturally from our above discussion. We set $`\mathrm{\Delta }=(a_2/a_4)^{1/2}\overline{\mathrm{\Delta }}`$, $`𝐐=(a_2/a_4)^{1/2}\overline{𝐪}`$ , $`a_0=A_0a_2^2/a_4`$ and $`\mathrm{\Omega }=(a_2^3/a_4^2)F`$. This leads to rewrite Eq.(7) for the free energy as :

$`F={\displaystyle \underset{𝐪}{}}|\overline{\mathrm{\Delta }}_𝐪|^2[A_0{\displaystyle \frac{1}{3}}\overline{q}^2+{\displaystyle \frac{1}{5}}\overline{q}^4]{\displaystyle \underset{𝐪_i}{}}\overline{\mathrm{\Delta }}_{𝐪_1}\overline{\mathrm{\Delta }}_{𝐪_2}^{}\overline{\mathrm{\Delta }}_{𝐪_3}\overline{\mathrm{\Delta }}_{𝐪_4}^{}[{\displaystyle \frac{1}{4}}{\displaystyle \frac{1}{12}}(\overline{q}_1^2+\overline{q}_2^2+\overline{q}_3^2+\overline{q}_4^2+\overline{𝐪}_1.\overline{𝐪}_3+\overline{𝐪}_2.\overline{𝐪}_4)]`$ (13)
$`+{\displaystyle \frac{1}{8}}{\displaystyle \underset{𝐪_i}{}}\overline{\mathrm{\Delta }}_{𝐪_1}\overline{\mathrm{\Delta }}_{𝐪_2}^{}\overline{\mathrm{\Delta }}_{𝐪_3}\overline{\mathrm{\Delta }}_{𝐪_4}^{}\overline{\mathrm{\Delta }}_{𝐪_5}\overline{\mathrm{\Delta }}_{𝐪_6}^{}`$ (14)
It is clear from this rescaling transformation and the resulting expression that, in the vicinity of the TCP (where $`a_2`$ is small), it is unnecessary to go beyond our sixth order expansion in $`\mathrm{\Delta }`$ and $`𝐪`$. It is also of interest to rewrite this free energy as a functional of $`\mathrm{\Delta }(𝐫)`$ by Fourier transform. This gives, after by parts integrations :

$`F={\displaystyle 𝑑𝐫[A_0|\overline{\mathrm{\Delta }}|^2\frac{1}{3}|\overline{\mathrm{\Delta }}|^2+\frac{1}{5}|^2\overline{\mathrm{\Delta }}|^2]}{\displaystyle 𝑑𝐫[\frac{1}{4}|\overline{\mathrm{\Delta }}|^4\frac{1}{24}[2(|\overline{\mathrm{\Delta }}|^2)^2+3(\overline{\mathrm{\Delta }}^2)(\overline{\mathrm{\Delta }}^2)]]}+{\displaystyle \frac{1}{8}}{\displaystyle 𝑑𝐫|\overline{\mathrm{\Delta }}|^6}`$ (15)
## IV THE LO SUBSPACE

As emphasized by LO all the states corresponding to the same wavector $`\overline{q}_0`$ but with different orientation for $`\overline{𝐪}_0`$ are degenerate right on the FFLO transition line. This degeneracy is lifted, at least partially, when one goes into the superfluid phase because of the coupling between the various plane waves produced by the nonlinear terms in the free energy. When one investigates the states selected in this process, one has to consider the subspace:

$`\overline{\mathrm{\Delta }}(𝐫)={\displaystyle }\mathrm{\Delta }_{𝐪_0}\mathrm{exp}(i\overline{𝐪}_0.𝐫)`$ (16)
of all the order parameters generated by these plane waves. We call this the LO subspace. In this section we will restrict to this subspace our search for the state appearing at the transition : we will look for the minimum of the free energy within this LO subspace. Actually LO looked for a lattice as a solution and restricted themselves to this kind of order parameter. However there is physically no basic reason to enforce this type of restriction. One could look for incommensurate structures or quasicrystal-like solutions. Even if these are not the lowest energy solution, they might be of interest as local solutions corresponding physically to defects. Therefore we have not set a periodicity condition on the solutions we have considered. Nevertheless let us indicate at once that the energetically favored solutions we have found within the LO subspace are actually periodic. Although considering only the LO subspace is an important restriction, it does not make the problem easy at all, although we solve it completely below.

Let us first show that, within this subspace, the standard FFLO transition line is not stable, and the transition is actually first order. With our reduced units the LO subspace corresponds to $`\overline{𝐪}_0^2=\frac{5}{6}`$. This minimizes the coefficient of the second order term in Eq.(9). The FFLO transition line is then given by $`A_0=\frac{5}{36}`$, which makes the second order term zero within the LO subspace. Let us then look at the fourth order term in Eq.(9). This amounts to calculate the functional derivative of the free energy on the FFLO transition line. In a standard second order phase transition it should always be positive, forcing the order parameter to be zero on the transition line. However in the present case it is not obvious that this is systematically so because of the interplay of the wavevectors in this term. Specifically we introduce a parameter $`\beta `$ which describes this effect for any fixed order parameter $`\mathrm{\Delta }(𝐫)`$. It is defined by :

$`2\beta \overline{q}_0^2{\displaystyle \underset{𝐪_i}{}}\overline{\mathrm{\Delta }}_{𝐪_1}\overline{\mathrm{\Delta }}_{𝐪_2}^{}\overline{\mathrm{\Delta }}_{𝐪_3}\overline{\mathrm{\Delta }}_{𝐪_4}^{}=\overline{q}_0^2{\displaystyle \underset{𝐪_i}{}}(\widehat{𝐪}_1.\widehat{𝐪}_3+\widehat{𝐪}_2.\widehat{𝐪}_4)\overline{\mathrm{\Delta }}_{𝐪_1}\overline{\mathrm{\Delta }}_{𝐪_2}^{}\overline{\mathrm{\Delta }}_{𝐪_3}\overline{\mathrm{\Delta }}_{𝐪_4}^{}={\displaystyle }\overline{\mathrm{\Delta }}^2(\overline{\mathrm{\Delta }}^{})^2+c.c.`$ (17)
For the simple case of the FF solution we have merely $`\beta =1`$. However this is the highest possible value and we can think of decreasing it, or even making it negative, by a proper choice of the order parameter in the LO space, although naturally $`\beta 1`$. Then the fourth order term in Eq.(9) is given by $`[\frac{1}{6}(\beta +2)\overline{q}_0^2\frac{1}{4}]|\overline{\mathrm{\Delta }}|^4`$ and has thus the same sign as $`5\beta +1`$ when we take $`\overline{q}_0^2=\frac{5}{6}`$. Hence for any order parameter with $`\beta \frac{1}{5}`$ the fourth order term is negative. We will find actually many such states. Now when we are on the FFLO transition line (the second order term is zero) and the fourth order term is negative, we decrease the free energy and make it negative just by taking the order parameter to be small and nonzero, which shows that the standard FFLO transition line is unstable. Naturally the larger the order parameter, the lower the free energy, but we have to stay in the range where the second and fourth order terms are the only ones important in our expansion, which means that the sixth order term is negligible. Then since we have a negative free energy, we can raise it back to zero by increasing $`A_0`$ beyond its value $`\frac{5}{36}`$ on the FFLO line, which means we go beyond this line in the $`(\overline{\mu },T)`$ phase diagram. In this way the second order term becomes positive. Eventually we will be limited by the growth of the sixth order term. The transition we have found corresponds to positive second and sixth order terms and a negative fourth order one. The free energy becomes negative for a nonzero value of the order parameter. We have thus found a first order transition beyond the standard FFLO transition line. This discussion about the order of the transition is the exact analogue of the one we made in the introduction for the $`T=0`$ situation.

Naturally it is of interest to minimize $`\beta `$ since it is rather natural to expect that the states corresponding to the minimum will lead to the stronger instability toward the first order transition. We show now that $`\beta \frac{1}{3}`$, the equality $`\beta =\frac{1}{3}`$ being obtained for any real order parameter. We make use of :

$`{\displaystyle }[\overline{\mathrm{\Delta }}^2(\overline{\mathrm{\Delta }}^{})^2|\overline{\mathrm{\Delta }}|^2|\overline{\mathrm{\Delta }}|^2]+c.c.={\displaystyle }[\overline{\mathrm{\Delta }}\overline{\mathrm{\Delta }}^{}c.c.]^20`$ (18)
where the equality occurs only for a real order parameter (within an irrelevant overall constant phase factor), and :

$`{\displaystyle }[\overline{\mathrm{\Delta }}^2(\overline{\mathrm{\Delta }}^{})^2+2|\overline{\mathrm{\Delta }}|^2|\overline{\mathrm{\Delta }}|^2]+c.c.={\displaystyle }|\overline{\mathrm{\Delta }}|^2(\overline{\mathrm{\Delta }}^2\overline{\mathrm{\Delta }}^{}+c.c.)=2\overline{q}_0^2{\displaystyle }|\overline{\mathrm{\Delta }}|^4`$ (19)
where the last step makes specific use of the form Eq.(11) for the order parameter. We have then :

$`2\beta \overline{q}_0^2{\displaystyle }|\overline{\mathrm{\Delta }}|^4={\displaystyle }\overline{\mathrm{\Delta }}^2(\overline{\mathrm{\Delta }}^{})^2+c.c.{\displaystyle }[\overline{\mathrm{\Delta }}^2(\overline{\mathrm{\Delta }}^{})^2{\displaystyle \frac{2}{3}}[\overline{\mathrm{\Delta }}^2(\overline{\mathrm{\Delta }}^{})^2|\overline{\mathrm{\Delta }}|^2|\overline{\mathrm{\Delta }}|^2]]+c.c.`$ (20)
$`={\displaystyle \frac{1}{3}}{\displaystyle }[\overline{\mathrm{\Delta }}^2(\overline{\mathrm{\Delta }}^{})^2+2|\overline{\mathrm{\Delta }}|^2|\overline{\mathrm{\Delta }}|^2]+c.c.={\displaystyle \frac{2}{3}}\overline{q}_0^2{\displaystyle }|\overline{\mathrm{\Delta }}|^4`$ (21)
hence $`\beta \frac{1}{3}`$.

Together with $`\beta `$ it is also intuitively convenient to consider $`\gamma `$ defined by :

$`\gamma {\displaystyle 𝑑𝐫|\overline{\mathrm{\Delta }}|^4}=\overline{q}_0^2{\displaystyle 𝑑𝐫[|\overline{\mathrm{\Delta }}|^2]^2}={\displaystyle \underset{𝐪_i}{}}(\widehat{𝐪}_1\widehat{𝐪}_2)^2\overline{\mathrm{\Delta }}_{𝐪_1}\overline{\mathrm{\Delta }}_{𝐪_2}^{}\overline{\mathrm{\Delta }}_{𝐪_3}\overline{\mathrm{\Delta }}_{𝐪_4}^{}`$ (22)
which is easily seen to satisfy $`\gamma =1\beta `$ from Eq.(14). We are thus interested in maximizing $`\gamma `$. From the first expression in Eq.(16) it is intuitively clear that $`\gamma `$ will be large when the (properly normalized) order parameter has strong spatial variations. In particular it will be better for $`|\overline{\mathrm{\Delta }}(𝐫)|^2`$ to have many nodes. This is more easily achieved if $`\overline{\mathrm{\Delta }}(𝐫)`$ has no imaginary part since in this case one has only to require that the real part is zero. This makes intuitively reasonable that $`\gamma `$ is maximized by a real order parameter. One can also come to this conclusion from the Fourier expansion in Eq.(16) (or also from Eq.(12)) : it is of interest to have as often as possible opposite wavevectors so that $`(\widehat{𝐪}_1\widehat{𝐪}_2)^2`$ takes as much as possible its maximum value, namely 4 . One is thus naturally led to an order parameter which is a combination of $`\mathrm{cos}(𝐪_0.𝐫+\phi _{𝐪_0})`$ with real coefficients, the directions of the $`𝐪_0`$ ’s being free. This corresponds merely to require that $`\overline{\mathrm{\Delta }}(𝐫)`$ is real in Eq.(11) by taking $`\mathrm{\Delta }_{𝐪_0}=\mathrm{\Delta }_{𝐪_0}^{}`$. On the other hand, since any real order parameter gives the maximum $`\gamma `$, it is not necessary that these cosines have equal weight.

We can now come back to our free energy Eq.(9) and find the best solution within the LO subspace, since we have found that the minimum $`\beta `$ is $`1/3`$, as soon as the order parameter is real, which implies $`\mathrm{\Delta }_{𝐪_0}=\mathrm{\Delta }_{𝐪_0}^{}`$ in Eq.(11). We introduce a measure $`\overline{\mathrm{\Delta }}`$ of the amplitude of the order parameter by setting :

$`{\displaystyle 𝑑𝐫|\overline{\mathrm{\Delta }}|^2}={\displaystyle \underset{𝐪_i}{}}\overline{\mathrm{\Delta }}_{𝐪_1}\overline{\mathrm{\Delta }}_{𝐪_2}^{}=N_2\overline{\mathrm{\Delta }}^2`$ (23)
where by definition $`N_2N`$ is the number of plane waves coming in Eq.(11). If all the planed waves have same amplitude, $`\overline{\mathrm{\Delta }}`$ is just the common value of these amplitudes. Then we define $`N_4`$ and $`N_6`$ by :

$`{\displaystyle |\overline{\mathrm{\Delta }}|^4}={\displaystyle \underset{𝐪_i}{}}\overline{\mathrm{\Delta }}_{𝐪_1}\overline{\mathrm{\Delta }}_{𝐪_2}^{}\overline{\mathrm{\Delta }}_{𝐪_3}\overline{\mathrm{\Delta }}_{𝐪_4}^{}=N_4\overline{\mathrm{\Delta }}^4`$ (24)
and :

$`{\displaystyle |\overline{\mathrm{\Delta }}|^6}={\displaystyle \underset{𝐪_i}{}}\overline{\mathrm{\Delta }}_{𝐪_1}\overline{\mathrm{\Delta }}_{𝐪_2}^{}\overline{\mathrm{\Delta }}_{𝐪_3}\overline{\mathrm{\Delta }}_{𝐪_4}^{}\overline{\mathrm{\Delta }}_{𝐪_5}\overline{\mathrm{\Delta }}_{𝐪_6}^{}=N_6\overline{\mathrm{\Delta }}^6`$ (25)
For the simple plane wave considered in Eq.(8), we had $`N_2=N_4=N_6=1`$. For the real order parameter we are interested in, the set of wavevectors $`\{𝐪_i\}`$ is made of N/2 pairs. If all the planed waves have same amplitude, one finds  for example by simple counting that $`N_4=3N(N1)`$ and $`N_6=5N(3N^29N+8)`$. Actually, once it is recognized that, from the counting procedure, $`N_4`$ and $`N_6`$ are polynomials in $`N`$ of order 2 and 3 respectively, the coefficient may easily be found by considering the cases $`N=2,4,6`$. With these notations the free energy Eq.(9) reduces to :

$`F=N_2\overline{\mathrm{\Delta }}^2[A_0{\displaystyle \frac{1}{3}}\overline{q}^2+{\displaystyle \frac{1}{5}}\overline{q}^4]N_4\overline{\mathrm{\Delta }}^4[{\displaystyle \frac{1}{4}}{\displaystyle \frac{\alpha }{2}}\overline{q}^2]+{\displaystyle \frac{1}{8}}N_6\overline{\mathrm{\Delta }}^6`$ (26)
where we have set $`\alpha =\frac{\beta +2}{3}`$ (for the simple plane wave considered in Eq.(8), we had $`\alpha =1`$).

We proceed now as we have done for Eq.(8). Minimizing F with respect to $`\overline{q}^2`$ we find for the extremum the condition $`\overline{q}^2=\frac{5}{6}\frac{5}{4}\alpha (N_4/N_2)\overline{\mathrm{\Delta }}^2`$, which implies that, in our considerations, we have for $`\overline{\mathrm{\Delta }}^2`$ an upper bound $`\overline{\mathrm{\Delta }}_{max}^2=(2/3\alpha )N_2/N_4`$. This leads, for the value of the free energy F at this extremum, to :

$`{\displaystyle \frac{F}{\overline{\mathrm{\Delta }}^2}}=N_2(A_0{\displaystyle \frac{5}{36}})+N_4\overline{\mathrm{\Delta }}^2({\displaystyle \frac{5\alpha }{12}}{\displaystyle \frac{1}{4}})+\overline{\mathrm{\Delta }}^4({\displaystyle \frac{N_6}{8}}{\displaystyle \frac{5\alpha ^2}{16}}{\displaystyle \frac{N_4^2}{N_2}})`$ (27)
Now the standard FFLO solution corresponds to $`A_0=\frac{5}{36}`$. Since we are interested in a better solution we want $`A_0>\frac{5}{36}`$ which makes the first term of $`F/\overline{\mathrm{\Delta }}^2`$ positive. On the other hand, for $`\overline{\mathrm{\Delta }}=\overline{\mathrm{\Delta }}_{max}`$, the sum of the last two terms in the right hand side (r.h.s.) of Eq.(21) can be written as $`\overline{\mathrm{\Delta }}_{max}^2(N_4/24\alpha )[5(\alpha \frac{3}{5})^2+2N_2N_6/N_4^2\frac{9}{5}]`$. This is always positive since we have $`N_2N_6/N_4^21`$ ( this results directly from $`𝑑𝐫|\overline{\mathrm{\Delta }}|^2.𝑑𝐫|\overline{\mathrm{\Delta }}|^6[𝑑𝐫|\overline{\mathrm{\Delta }}|^4]^2`$). If we assume that the second term in the r.h.s. of Eq.(21) is positive, this implies that the free energy is always positive. Therefore if we want to find non positive values for F it is necessary to have a negative coefficient for the second term in the r.h.s. of Eq.(21), which means $`\alpha <\frac{3}{5}`$. This is possible since the minimum $`\beta =\frac{1}{3}`$ we have found above corresponds to a minimum $`\alpha =\frac{5}{9}`$.

Then the quadratic form Eq.(21) can be zero if we meet the condition $`4(A_0\frac{5}{36})(2\frac{N_2N_6}{N_4^2}5\alpha ^2)(\frac{5\alpha }{3}1)^2`$. This leads to the following result for the transition line to the FFLO phase :

$`A_0={\displaystyle \frac{5}{36}}+{\displaystyle \frac{1}{8}}{\displaystyle \frac{(1\frac{5\alpha }{3})^2}{\frac{N_2N_6}{N_4^2}\frac{5}{2}\alpha ^2}}`$ (28)
In the $`(\overline{\mu }T)`$ plane this line is higher than the standard second order FFLO transition line, which is given by $`A_0=\frac{5}{36}`$. On the other hand the value of $`\overline{\mathrm{\Delta }}_m`$ which gives $`F=0`$ at the threshold given by Eq.(21) is :

$`\overline{\mathrm{\Delta }}_m^2={\displaystyle \frac{N_2}{N_4}}{\displaystyle \frac{1\frac{5\alpha }{3}}{\frac{N_2N_6}{N_4^2}\frac{5}{2}\alpha ^2}}`$ (29)
which is clearly nonzero and our transition is quite explicitely a first order transition (note that the condition $`\overline{\mathrm{\Delta }}_m<\overline{\mathrm{\Delta }}_{max}`$ is necessarily satisfied since F is zero for $`\overline{\mathrm{\Delta }}_m`$ and positive for $`\overline{\mathrm{\Delta }}_{max}`$ ).

Coming back to the location of the transition line Eq.(23) we consider now how to optimize it. First we see from the numerator that it is advantageous to lower $`\alpha `$ as much as possible. Indeed our minimum $`\alpha =\frac{5}{9}`$ is quite close to the limiting value $`\alpha =\frac{3}{5}`$ so the variation of the denominator with $`\alpha `$ is irrelevant. Hence we are lead to make the fourth order term in the free energy Eq.(9) as negative as possible by minimizing $`\beta `$, as we anticipated at the beginning of section IV. Once we have taken $`\alpha =\frac{5}{9}`$, we see that it is of interest to take $`N_2N_6/N_4^2`$ as small as possible (we have seen that it is bounded from below by $`1`$). From their definitions Eq.(18) and Eq.(19) we can evaluate $`N_4`$ and $`N_6`$ for a general order parameter Eq.(11) $`\overline{\mathrm{\Delta }}(𝐫)=\mathrm{\Delta }_{𝐪_i}\mathrm{exp}(i𝐪_i.𝐫)=2|\mathrm{\Delta }_{𝐪_i}|\mathrm{cos}(𝐪_i.𝐫+\phi _{𝐪_i})`$. One finds :

$`N_4\overline{\mathrm{\Delta }}^4=3({\displaystyle \underset{i}{}}|\mathrm{\Delta }_{𝐪_i}|^2)^23{\displaystyle \underset{i}{}}|\mathrm{\Delta }_{𝐪_i}|^4`$ (30)
and :

$`N_6\overline{\mathrm{\Delta }}^6=15({\displaystyle \underset{i}{}}|\mathrm{\Delta }_{𝐪_i}|^2)^345({\displaystyle \underset{i}{}}|\mathrm{\Delta }_{𝐪_i}|^2)^2{\displaystyle \underset{i}{}}|\mathrm{\Delta }_{𝐪_i}|^4+40{\displaystyle \underset{i}{}}|\mathrm{\Delta }_{𝐪_i}|^6`$ (31)
Actually one can recognize that, for symmetry reasons, the results involve only $`_i|\mathrm{\Delta }_{𝐪_i}|^n`$ with $`n=2,4,6`$, because odd powers of cosines average to zero. Hence the results assume necessarily the general form given by Eq.(24) and Eq.(25). Then the coefficients are easily obtained from the specific case where all the amplitudes $`|\mathrm{\Delta }_{𝐪_i}|`$ are equal.

Since the wavevectors are paired it is now more convenient to sum over pairs from now on. Defining $`a_i=|\mathrm{\Delta }_{𝐪_i}|^2/_i|\mathrm{\Delta }_{𝐪_i}|^2`$ (implying $`_ia_i=1`$) we have :

$`{\displaystyle \frac{9}{10}}N_2N_6/N_4^2={\displaystyle \frac{69S_2+4S_3}{(2S_2)^2}}`$ (32)
where $`S_2=_ia_i^2`$ and $`S_3=_ia_i^3`$. When we have a single pair $`S_2=S_3=1`$ and the r.h.s. of Eq.(26) is equal to 1. We show now that it is otherwise larger than 1. Since $`(2S_2)^2<43S_2`$ because $`S_2<1`$ when we have more than a single pair, it is enough to prove that $`3(1S_2)2(1S_3)`$. This is in turn verified because we can write for the left-hand side $`1S_2=(_ia_i)^2_ia_i^2=2_{i<j}a_ia_j`$. In a similar way we have in the right-hand side $`1S_3=(_ia_i)^3_ia_i^3=3_{i<j}a_ia_j(a_i+a_j)+6_{i<j<k}a_ia_ja_k=3_{i<j}a_ia_j9_{i<j<k}a_ia_ja_k`$. So our statement is correct since it is just equivalent to $`_{i<j<k}a_ia_ja_k0`$ (the equality holds when we have only two pairs since one can not have three different indices). Therefore we come to the conclusion that $`N_2N_6/N_4^2`$ is minimized when we take a single pair of plane waves with wavevectors ($`𝐪,𝐪`$), corresponding to a simple order parameter proportional to $`\mathrm{cos}(𝐪.𝐫)`$, which has hence a planar symmetry. It is unfavourable to increase the number of plane waves. This is more easily seen in the particular situation where these $`N`$ plane waves have the same amplitude. In this case we merely have $`N_2N_6/N_4^2=[15N^245N+40]/9(N1)^2`$ which increases regularly with increasing N, and is minimum for $`N=2`$ . In this last case $`N_2N_6/N_4^2=\frac{10}{9}`$. It is interesting to remark that this conclusion is opposite to what one would obtain by considering the fourth order term alone in Eq.(9) and omitting the sixth order one. Since the fourth order term grows (compared to the second order one) with the number $`N`$ plane waves, one would conclude that it is better to increase this number. The opposite turns out to be true because the sixth order term grows even faster with the number of plane waves. This shows quite clearly that, in contrast with what one might hope, the consideration of the fourth order term is not enough to conclude about the actual ground state of system.

We find then explicitely $`A_0=\frac{5}{36}+\mathrm{2.02\hspace{0.25em}10}^3`$. We note that this is quite close to the standard FFLO transition itself. For comparison the standard Clogston-Chandrasekhar  first order transition is given, as we have seen, by $`A_0=\frac{1}{8}`$ so the difference with the standard FFLO transition is $`\mathrm{1.39\hspace{0.25em}10}^2`$ which is also rather small. Since we know that the standard Clogston-Chandrasekhar transition and standard FFLO transition stay close beyond the vicinity of the TCP and that this proximity between these two lines extends down to zero temperature, it seems quite possible that the same is true for our first order transition line. This is indeed what has been found numerically by Matsuo et al. . On the other hand we know from the work of Burkhardt and Rainer  that the transition is second order in a 2D situation. So it is of some interest to consider formally an arbitrary dimension to see how one goes from the 2D to the 3D situation. For the case of our planar order parameter, there is no difficulty in working with an arbitrary dimension D. As we have mentionned at the end of section II, one has just to modify some coefficients in Eq.(7) : here one has to multiply the $`Q^2`$ terms by $`3/D`$ and the $`Q^4`$ terms by $`15/D(D+2)`$. Reproducing then the above analysis, we obtain that the FFLO transition occurs for $`A_0=\frac{D+2}{12D}`$ with a wavevector $`\overline{q}^2=\frac{D+2}{6}`$. One finds again that, in order to obtain a phase transition higher than FFLO, the coefficient of the second term in the r.h.s. of Eq.(21) has to be negative. But this coefficient is now in general $`N_4(\frac{\alpha (D+2)}{D}1)/4`$ with explicitely $`\alpha =\frac{5}{9}`$. So we find $`D_c=2.5`$ as the critical dimension to have a first order transition, which is halfway between the physical situations $`D=2`$ and $`D=3`$. Hence, in contrast to what one might think, the first order transition we have found is not ’very near to be second order’, despite the proximity between the first order transition line and the standard second order FFLO line. The result Eq.(21) for the location of the transition for this planar state becomes :

$`A_0={\displaystyle \frac{D+2}{12D}}+{\displaystyle \frac{(2D5)^2}{15D(7D10)}}`$ (33)
with the value of $`\overline{\mathrm{\Delta }}_m`$ of $`\overline{\mathrm{\Delta }}`$ at the threshold given by :

$`\overline{\mathrm{\Delta }}_m^2={\displaystyle \frac{4}{5}}{\displaystyle \frac{2D5}{7D10}}`$ (34)
In two dimensions we can again make use of Eq.(20) to compare all the possible states. As we have seen the transition is second order and we can neglect the sixth order term as in the original LO analysis. The fourth order term is positive and, roughly speaking, below the transition the state with smaller fourth order term is selected. This leads to take $`N_4/N_2^2=3(2S_2)/2`$ as small as possible. So we find again that the single pair state $`N=2`$ is selected, which corresponds to the state investigated by Burkhardt and Rainer .

In closing this section, let us note that we have have found a non- degenerate minimum for the free energy, namely the planar state $`\mathrm{cos}(𝐪.𝐫)`$ corresponding to two plane waves with opposite wavevectors and equal amplitude in Eq.(11). However by adding to this order parameter another cosine with a small amplitude, namely $`ϵ\mathrm{cos}(𝐪_\mathrm{𝟏}.𝐫)`$, we can have a free energy arbitrarily close to the planar state free energy. Conversely the location for the transition can be made arbitrarily close to the planar state transition. These states are highly degenerate since the relative directions of $`𝐪`$ and $`𝐪_\mathrm{𝟏}`$ do not matter. These are some kind of excited states for our system with arbitrarily small energy and it is possible to speculate that they play a role in the physics of our system.

## V BEYOND THE LO SUBSPACE

In the preceding section we have restricted to the LO subspace our search for the state with lowest free energy at the transition and therefore with highest critical temperature. While this restriction is justified when the free energy expansion is limited to the fourth order term (this term being treated as a perturbation), it is no longer valid when the sixth order term is included. And indeed the actual minimum state does not belong to the LO subspace, since it can be checked that Eq.(11) can not satisfy the Euler-Lagrange equation corresponding to Eq.(10) because of the non linear terms. However we have found in the preceding section a non degenerate minimum corresponding to the planar state. The selection of the wavevectors of this state was due to the fourth order term while the number of components was controled by the sixth order one, which is otherwise quite structureless. If the nonlinearities produce only a rather small change on the solution, it is reasonable to look for the actual solution as being ’near’ the solution $`\mathrm{cos}(𝐪.𝐫)`$ we have found in section IV. In particular it is reasonable to look for an order parameter with a one-dimensional dependence on $`𝐫`$ and we will also assume it to be real. And indeed we will find a solution which is quite near $`\mathrm{cos}(𝐪.𝐫)`$. Although this argumentation is only a self-consistent one (we can not exclude that a strong nonlinear modification of a three-dimensional solution produces the actual minimum), we note that the numerical exploration of Houzet et al  has indeed produced a real one-dimensional order parameter for the minimum.

With a one-dimensional real order parameter the reduced free energy simplifies into :

$`F={\displaystyle 𝑑x[A_0\overline{\mathrm{\Delta }}^2\frac{1}{3}\overline{\mathrm{\Delta }}^{}_{}{}^{}2+\frac{1}{5}\overline{\mathrm{\Delta }}^{{}_{}{}^{\prime \prime }2}]}{\displaystyle 𝑑x[\frac{1}{4}\overline{\mathrm{\Delta }}^4\frac{5}{6}\overline{\mathrm{\Delta }}^2\overline{\mathrm{\Delta }}^{}_{}{}^{}2]}+{\displaystyle \frac{1}{8}}{\displaystyle 𝑑x\overline{\mathrm{\Delta }}^6}`$ (35)
where $`\overline{\mathrm{\Delta }}^{}`$ and $`\overline{\mathrm{\Delta }}^{\prime \prime }`$ are first and second derivative of $`\overline{\mathrm{\Delta }}`$ with respect to $`x`$. Although we deal with a nonlinear problem we can still minimize with respect to the amplitude of the order parameter, just as we have done at the end of the preceding section (this works actually also for a three-dimensional order parameter). We set $`\overline{\mathrm{\Delta }}(x)=a\delta (x)`$, where $`\delta (x)`$ is a normalized spatial function (for example by $`\delta ^2=1`$ ). $`\delta (x)`$ gives the shape of the order parameter while $`a`$ corresponds to its amplitude. When we substitute in Eq.(29) we obtain for $`F/a^2`$ a quadratic form in $`a^2`$. Writing again that, for a specific function $`\delta (x)`$, this form has double root at the transition leads us to an expression for $`A_0`$ which does not depend anymore on the normalization of $`\delta `$ :

$`A_0={\displaystyle \frac{[\delta ^4\frac{10}{3}(\delta \delta ^{})^2]^2}{8\delta ^2\delta ^6}}+{\displaystyle \frac{\frac{1}{3}\delta ^{}_{}{}^{}2\frac{1}{5}\delta ^{{}_{}{}^{\prime \prime }2}}{\delta ^2}}`$ (36)
and we want now to maximize A<sub>0</sub>. Similarly the condition that the free energy is zero at the transition imply :

$`{\displaystyle \overline{\mathrm{\Delta }}^6}={\displaystyle \overline{\mathrm{\Delta }}^4}{\displaystyle \frac{10}{3}}\overline{\mathrm{\Delta }}^2\overline{\mathrm{\Delta }}^{}_{}{}^{}2`$ (37)
or equivalently :

$`A_0{\displaystyle \overline{\mathrm{\Delta }}^2}={\displaystyle \frac{1}{8}\overline{\mathrm{\Delta }}^6}+{\displaystyle \frac{1}{3}}\overline{\mathrm{\Delta }}^{}_{}{}^{}2{\displaystyle \frac{1}{5}}\overline{\mathrm{\Delta }}^{{}_{}{}^{\prime \prime }2}`$ (38)
which gives the amplitude of the order parameter.

The condition that $`A_0`$ is maximum can be obtained from Eq.(30) but it is easier to deduce it from Eq.(29). One finds the ordinary nonlinear differential equation :

$`{\displaystyle \frac{1}{5}}\overline{\mathrm{\Delta }}^{\prime \prime \prime \prime }+({\displaystyle \frac{1}{3}}{\displaystyle \frac{5}{6}}\overline{\mathrm{\Delta }}^2)\overline{\mathrm{\Delta }}^{\prime \prime }{\displaystyle \frac{5}{6}}\overline{\mathrm{\Delta }}\overline{\mathrm{\Delta }}^{}_{}{}^{}2+{\displaystyle \frac{3}{8}}\overline{\mathrm{\Delta }}^5{\displaystyle \frac{1}{2}}\overline{\mathrm{\Delta }}^3+A_0\overline{\mathrm{\Delta }}=0`$ (39)
(analytically this equation can be integrated once and then reduced to a second order nonlinear differential equation but this is of no real help for numerics). One checks readily that the cosine form of the order parameter $`\overline{\mathrm{\Delta }}(x)=a\mathrm{cos}(\overline{q}x)`$ is not a solution. It satisfies only the linear part of the equation, which leads to A$`{}_{0}{}^{}=\overline{q}^2/3\overline{q}^4/5`$. Looking for the maximum of $`A_0`$ gives the FFLO result $`A_0=5/360.1389`$ with $`\overline{q}_{FFLO}^2=5/6`$. This FFLO solution corresponds to keep only the second term in the r.h.s. of Eq.(30). This term produces for $`A_0`$ a kind of (inverted) effective potential $`\overline{q}^2/3\overline{q}^4/5`$ which has a strong maximum for the FFLO solution. When we write a Fourier expansion of the solution this potential selects quite effectively the wavevectors in the close vicinity of the FFLO result $`\overline{q}_{FFLO}`$.

Indeed the numerical exploration of Eq.(33) for $`A_0`$ very close to the maximum gives for the solutions one or two basis frequencies $`q_0`$ and $`q_1`$ which are close to $`\overline{q}_{FFLO}`$. The other frequencies appearing in a Fourier analysis are simply odd combinations of $`q_0`$ and $`q_1`$ like $`2q_0\pm q_1`$, $`2q_1\pm q_0`$, $`3q_0`$, $`4q_1\pm q_0`$, etc. The weights of these frequencies depend on their order, i.e. the higher the frequency, the smaller the weight. This shows explicitely the strength of the effective FFLO potential in A<sub>0</sub> (the second term in the r.h.s. of Eq.(30)). We have investigated analytically the efficiency of a small frequency splitting by writing $`q_1=q_0+ϵ`$ with $`ϵ0`$ (but $`ϵ0`$ because this limit is singular). However a single frequency $`q_1=q_0`$ turned out to be always better. We are finally lead to the conclusion that going beyond the LO subspace produces only small corrections due to nonlinearities. These corrections correspond to odd harmonics of the fundamental frequency $`\overline{q}`$ and they are small because of the efficiency of the effective FFLO potential. Actually if we assume from the start that the corrections to $`\delta _0(x)=\mathrm{cos}(\overline{q}x)`$ are small, we can write $`\delta (x)=\delta _0(x)+\delta _1(x)`$, with $`\delta _1(x)`$ small, and perform a first order expansion of Eq.(30). One finds readily that it is most favorable to take $`\delta _1(x)`$ proportional to $`\mathrm{cos}^3(\overline{q}x)`$ (another contribution from $`\mathrm{cos}^5(\overline{q}x)`$ is quite small), but one has to go to second order to find the amplitude. Actually, once it is proved that this harmonic expansion is the correct answer, it is much easier to minimize Eq.(30) numerically which avoids cumbersome calculations. Specifically we considered the trial function $`\delta (x)=\mathrm{cos}(\overline{q}x)+a_3\mathrm{cos}(3\overline{q}x+\varphi _1)+a_5\mathrm{cos}(5\overline{q}+\varphi _2)`$ and numerically maximized A<sub>0</sub> (some more complete trial functions that we also tried eventually reduced to this form when maximized). We found $`a_3=\mathrm{1.33.10}^2`$, $`a_5=\mathrm{1.62.10}^4`$, $`\overline{q}=0.793`$, $`\varphi _{1,2}=0`$ and $`A_0=0.141604`$. This form is very close to the cosine solution and the improvement in A<sub>0</sub> compared to our result Eq.(22) for the cosine solution $`A_0=31/2200.14091`$ is pretty small in absolute values, although it gives a significant increase of $`\mathrm{7\hspace{0.25em}10}^4`$ to our gain of $`\mathrm{2.02\hspace{0.25em}10}^3`$ compared to the FFLO result. Our result $`\overline{q}=0.793`$ has to be compared to our result in LO subspace from Eq.(23) $`\overline{q}=0.829`$ and to $`\overline{q}_{FFLO}=0.913`$.

## VI  CONCLUSION

In this paper we have explored analytically the nature of the transition to the FFLO superfluid phases in the vicinity of the tricritical point, where these phases begin to appear. This region is convenient for the analytical study we make because, in the vicinity of this point, one can make use of an expansion of the free energy up to sixth order, both in order parameter amplitude and in wavevector. Despite this simplification one has still a complex nonlinear problem to solve. We have first explored the minimization of this free energy within the LO subspace, made of arbitrary superpositions of plane waves. We have seen that the standard second order FFLO phase transition is unstable and that a first order transition occurs at higher temperature. Within this subspace we have shown that it is favorable to have a real order parameter. Then among these states we have shown that those with the smallest number of plane waves are prefered. This leads to retain only two plane waves, corresponding to an order parameter with a $`\mathrm{cos}(𝐪_0.𝐫)`$ dependence, in agreement with preceding work . Finally we have shown that, when releasing the constraint of working within the LO subspace, the order parameter at the transition is only very slightly modified by higher harmonics contributions and we have been able to ascribe this result to the strong selection of the wavevector caused by the second order terms of the free energy, corresponding physically to the standard FFLO transition.

We are very grateful to X. Leyronas for stimulating discussions and to K. Nagai for making us aware of his work, during a meeting very nicely organized by G. Eska and H. Brandt for D. Rainer’s 60th birthday.

* Laboratoire associé au Centre National de la Recherche Scientifique et aux Universités Paris 6 et Paris 7.