text:
# The Little Higgs from a Simple Group

## 1 Introduction

The Standard Model (SM) is well supported by all high energy data . Precision tests match predictions including one-loop quantum corrections. This suggest the SM is a valid description of Nature up to energies in the multi-TeV range with a Higgs mass which is less than about 200 GeV .

This picture is not satisfying because the Higgs gets one-loop quadratically divergent corrections to its squared mass. The most significant corrections come from loops of top quarks, $`W`$-bosons and the Higgs (Figure 1). The top loop is the most severe - demanding a contribution to the Higgs mass of 200 GeV or less requires a momentum cutoff $`\mathrm{\Lambda }_{top}\stackrel{<}{_{}}700`$ GeV. If one allows fine tuning of order 10% between this correction and a counter term, one still needs $`\mathrm{\Lambda }_{top}\stackrel{<}{_{}}2`$ TeV.

The need to cancel quadratic divergences along with the consistency of electroweak precision measurements with the SM suggests new weakly coupled physics at $``$ 1 TeV.

Supersymmetry softly broken around 1 TeV is an example of new physics that meets these criteria. Loops of superpartners cancel all quadratic divergences in the SM. The minimal supersymmetric standard model’s most compelling feature is its suggestive unification of couplings. Its least compelling feature is the fact that it has over 100 new parameters and yet the current bound on the Higgs mass requires fine-tuning of parameters for successful electroweak symmetry breaking. Models that improve this situation exist, and thus weakly coupled superpartners at the weak scale remains an interesting option. However, as we close in on the parameter space of the minimal supersymmetric standard model, it is of great interest to find alternative weakly coupled theories of electroweak symmetry breaking.

A new class of models, called “Little Higgs” theories \[3-8\], produce a light Higgs boson with weakly coupled physics up to 10s of TeV. As in composite Higgs models , the Higgs is a pseudo-Nambu-Goldstone boson (PNGB) and is massless at tree-level. Its mass is protected by a global symmetry which is spontaneously broken. The symmetry is explicitly broken by weakly coupled operators in the theory which become the normal SM couplings below 1 TeV. The Little Higgs trick is that no one operator alone explicitly breaks the global symmetry protecting the Higgs mass, and therefore no quadratically divergent contribution exists at one loop. The purpose of this article is to present a new model of this type with the simplest gauge group to date.

In the next section, we show that simply extending the electroweak gauge group to $`SU(3)\times U(1)`$ automatically removes one-loop quadratic divergences from the gauge and top loops. We’ll show that if two scalar triplets have vacuum expectation values (VEVs) in the same direction, breaking the gauge group to $`SU(2)\times U(1)`$, there exists an $`SU(2)`$ doublet whose mass is protected from quadratic divergences at one loop. Adding the minimum (weakly coupled) top sector to generate a Yukawa coupling automatically protects the Higgs mass at one loop. The key is that by having both triplets get VEVs, both produce a set of NGBs (including doublets), either of which could have been eaten by the massive gauge fields. It is the fact that no one operator in the tree-level Lagrangian contains both triplets which protects the doublet mass. The symmetry protecting the Higgs mass in this model is an approximate $`[SU(3)/SU(2)]^2`$. We discuss the non-linear sigma model version of this theory and show how to get there from the linear sigma model.

In Section 3 we present a complete model, including a quartic Higgs coupling, which has no quadratic divergences at one loop. Extending the gauge group to $`SU(4)\times U(1)`$ allows one to generate a quartic coupling through a vacuum misalignment mechanism. The non-linear sigma model is $`[SU(4)/SU(3)]^4`$. The result is a two Higgs doublet model with heavy $`SU(4)/SU(2)`$ gauge bosons and new scalars at a TeV.

In Section 4 we discuss the spectrum of the $`SU(4)`$ model and indicate some of the phenomenological constraints \[10-14\]. Because our model has no more gauge couplings than the SM, all couplings and masses in the gauge sector are a function of the breaking scale $`f`$ and measured gauge couplings. We also comment on flavor issues in these models.

In the last section we discuss possibilities for future directions.

## 2 Little Higgs from a simple gauge group

In this section we describe a very simple extension of the SM which keeps the Higgs naturally light by employing the little Higgs mechanism. The model contains new particles and couplings which cancel the quadratic divergences from the top quark loop and from the SM gauge interactions. This is accomplished by enlarging the $`SU(2)`$ weak gauge interactions to $`SU(3)`$. The cancellation of the divergences from the Higgs self-coupling is more difficult and is described in Section 4.

Our mechanism for eliminating the one-loop divergence from gauge interactions (Figure 1.b) can be understood as follows: an $`SU(3)`$ gauge group is spontaneously broken near 1 TeV by a vacuum expectation value (VEV) of two $`SU(3)`$ triplet scalars. The triplet expectation values are aligned so that both vevs leave the same $`SU(2)`$ unbroken. Ignoring their coupling trough the gauge interactions each scalar breaks a global $`SU(3)SU(2)`$, each yielding 5 Nambu-Goldstone bosons (NGBs). With the $`SU(3)`$ gauge interactions turned on, the diagonal linear combination of NGBs is eaten, but the orthogonal linear combination remains massless at tree level. The masslessness of these modes is easily understood by noticing that in absence of a direct coupling between the two scalar triplets each of them “thinks” that it is the only field which breaks the gauge group and therefore contains exact NGBs. Quantum corrections from loops involving the gauge bosons generate couplings between the two triplets and therefore a mass for the pseudo-Nambu-Goldstone bosons (PNGB). However, there is no quadratically divergent one-loop diagram involving both scalar triplets. Finite and log-divergent diagrams contribute small scalar masses of order $`g/4\pi f100`$ GeV.

More concretely, consider an $`SU(3)`$ gauge theory with two scalar fields transforming as (complex) triplets $`𝚽_i,i=1,2,`$ of the gauge group with a potential

$$\frac{\lambda ^2}{2}(𝚽_1^{}𝚽_1f^2)^2+\frac{\lambda ^2}{2}(𝚽_2^{}𝚽_2f^2)^2$$
(1)
which generates VEVs for the $`𝚽`$ s. For simplicity, we assume equal couplings and VEVs for both triplets. This defines a linear sigma model with two global $`SU(3)`$ symmetries acting on the two triplets. The spontaneous breaking $`[SU(3)]^2[SU(2)]^2`$ yields five NGBs from each scalar. We now weakly gauge an $`SU(3)`$ such that both scalars are triplets under the gauge symmetry. This explicitly breaks the $`[SU(3)]^2`$ global symmetry to diagonal $`SU(3)`$. After spontaneous symmetry breaking the five “diagonal” NGBs are eaten by the Higgs mechanism and the five orthogonal linear combinations are PNGBs. The symmetry which they correspond to – “axial” $`SU(3)`$ – is explicitly broken by the gauge couplings. But since the tree level scalar potential respects both $`SU(3)`$’s the PNGBs remain massless at tree level.

We parametrize the scalars as

$`𝚽_1=e^{i\mathrm{\Theta }_1/f}\left(\begin{array}{c}0\\ 0\\ f+\rho _1\end{array}\right)=e^{i\mathrm{\Theta }_{eaten}/f}e^{i\mathrm{\Theta }/f}\left(\begin{array}{c}0\\ 0\\ f+\rho _1\end{array}\right)`$ (8)
$`𝚽_2=e^{i\mathrm{\Theta }_2/f}\left(\begin{array}{c}0\\ 0\\ f+\rho _2\end{array}\right)=e^{i\mathrm{\Theta }_{eaten}/f}e^{i\mathrm{\Theta }/f}\left(\begin{array}{c}0\\ 0\\ f+\rho _2\end{array}\right)`$ (15)
where the first parametrization is the most obvious, but the second is more convenient because it separates the eaten modes $`\mathrm{\Theta }_{eaten}`$ from the PNGBs $`\mathrm{\Theta }`$. Note that $`\mathrm{\Theta }_{eaten}`$ shifts under diagonal (vector) $`SU(3)`$ transformations whereas $`\mathrm{\Theta }`$ shifts under “axial” $`SU(3)`$. The $`\rho _i`$ are radial modes which obtain masses $`m_\rho =\lambda f`$ from the potential. Since we are interested in the physics of the PNGBs we will suppress the eaten fields $`\mathrm{\Theta }_{eaten}`$. Furthermore we will take the limit $`\lambda 4\pi `$ in which the radial modes decouple and our linear sigma model turns into the corresponding non-linear sigma model. The non-linear sigma model obtained by integrating out the $`\rho _i`$ includes non-renormalizable interactions which become strongly coupled at $`\mathrm{\Lambda }=4\pi f`$; at this scale the non-linear sigma model description breaks down. For most of this paper we will use the non-linear sigma model description because it is more general. It focuses on the physics of the PNGBs. Details of the UV theory which lead to the $`[SU(3)]^2`$ symmetry breaking are encoded in higher dimensional operators and decouple from the relevant physics. The non-linear sigma model field are then parameterized as

$`\mathrm{\Phi }_1=e^{i\mathrm{\Theta }/f}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ f\hfill \end{array}\right),\mathrm{\Phi }_2=e^{i\mathrm{\Theta }/f}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ f\hfill \end{array}\right)`$ (22)
where we suppressed the eaten fields and removed the bold type to distinguish non-linear from linear sigma model fields.

Expanded out in components, the five PNGBs in $`\mathrm{\Theta }`$ are

$`\mathrm{\Theta }=\mathrm{\Theta }^aT^a={\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{cc}\begin{array}{cc}0\hfill & 0\hfill \\ 0\hfill & 0\hfill \end{array}& h\\ h^{}& 0\end{array}\right)+{\displaystyle \frac{\eta }{4}}\left(\begin{array}{ccc}\hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 1& \hfill 0\\ \hfill 0& \hfill 0& \hfill 2\end{array}\right).`$ (30)
The $`T^a`$ are the usual $`SU(3)`$ generators, $`a`$ runs from $`4\mathrm{}8`$, and normalizations were chosen such that $`h`$ and $`\eta `$ have canonical kinetic terms. Under the unbroken $`SU(2)`$ gauge symmetry $`h`$ transforms like the SM Higgs, i.e. it is a complex doublet, and $`\eta `$ is a neutral scalar.

Lets return to the linear sigma model. Since the gauge interactions explicitly break the axial $`SU(3)`$ symmetry which protects the PNGBs we expect that quantum corrections from gauge interactions will generate a mass for them. In the linear sigma model the $`SU(3)`$ gauge interactions are

$$\left|(_\mu +igA_\mu )𝚽_1\right|^2+\left|(_\mu +igA_\mu )𝚽_2\right|^2.$$
(31)
At one loop, there are two quadratically divergent diagrams which contribute to scalar masses (Figure 2). For the following arguments we choose $`_\mu A^\mu =0`$ gauge. This is convenient because diagrams involving the trilinear gauge couplings cannot contribute to the scalar potential (no derivatives) in this gauge. Thus the second diagram vanishes.

Diagrams contributing to the scalar potential are shown in Figure 3. The first diagram is quadratically divergent, the second one is log-divergent, and diagrams with even more $`𝚽`$ insertions would be finite.

The first diagram gives

$$\mathrm{}\frac{g^2\mathrm{\Lambda }^2}{16\pi ^2}(𝚽_1^{}𝚽_1+𝚽_2^{}𝚽_2),$$
(32)
which preserves both $`SU(3)`$ symmetries and renormalizes the $`𝚽_i`$ potentials, eq. (1). Again, taking $`\mathrm{\Lambda }=4\pi f`$ this contribution is not larger than the tree level terms already present in the potential and does not destabilize the desired vacuum. Thus at the level of quadratically divergent diagrams, the PNGBs remain massless.

Log-divergent and finite diagrams do contribute to PNGB masses. For example, the second diagram in Figure 3 generates

$`\mathrm{}{\displaystyle \frac{g^4}{16\pi ^2}}|𝚽_2^{}𝚽_1|^2\mathrm{log}(\mathrm{\Lambda }^2/f^2){\displaystyle \frac{f^2}{16\pi ^2}}h^{}h.`$ (33)
This operator explicitly violates axial $`SU(3)`$ and contains a mass for the PNGBs. However, the mass is of order $`f/4\pi `$ which is sufficiently small. Note that the mass-squared generated is positive, which means that this term stabilizes the vacuum with aligned expectation values for $`𝚽_1`$ and $`𝚽_2`$. The top quark loop which we discuss in the next section generates a negative mass-squared and thereby triggers electroweak symmetry breaking.

This analysis could have also been performed in the non-linear sigma model. The expression for the gauge couplings is identical to eq. (31) except that the $`𝚽_i`$ are replaced by $`\mathrm{\Phi }_i`$. Loops like Figure 3.a generate $`|\mathrm{\Phi }_i^{}\mathrm{\Phi }_i|=f^2`$, a quadratically divergent contribution to the cosmological constant but no mass. Finite and log-divergent diagrams generate $`|\mathrm{\Phi }_2^{}\mathrm{\Phi }_1|^2`$ which does contain a Higgs mass of order 100 GeV if $`f`$ 1 TeV.

The absence of quadratic divergences can also be understood by noting that the diagram in Figure 3.a only involves one of the non-linear sigma model fields. Thus it is identical to the corresponding diagram in a theory with only one $`\mathrm{\Phi }`$. But in this theory $`\mathrm{\Theta }`$ would be eaten by the Higgs mechanism, thus it cannot have any non-derivative couplings in the Lagrangian. Therefore, only diagrams which involve both $`\mathrm{\Phi }`$’s (Figure 3.b) can contribute to the Higgs mass but they necessarily involve more internal propagators and are not quadratically divergent.

### 2.1 Top Yukawa coupling

We now show that it is straightforward to add fermions and a top Yukawa coupling which does not upset the radiative stability of the Higgs mass. This can be done for both the linear and non-linear sigma models but we will only present the analysis of the non-linear model.

Since $`SU(2)`$-weak is embedded into an $`SU(3)`$ gauge group, the top-bottom $`SU(2)`$-doublet is enlarged to a triplet $`Q^T=(t,b,\chi )`$. A mass of order $`f`$ for the extra fermion $`\chi `$ in the triplet and Yukawa couplings for the top quark are generated from couplings to the $`\mathrm{\Phi }`$’s

$`_{top}=\lambda _1\chi _1^c\mathrm{\Phi }_1^{}Q+\lambda _2\chi _2^c\mathrm{\Phi }_2^{}Q,`$ (34)
where $`\chi _i^c`$ are Weyl fermions with the quantum numbers of the $`SU(2)`$-singlet component of the top quark.111Note that despite superficial similarities between eq. (34) and the top-color see-saw  there is an important difference. Unlike the top-$`\chi `$ Higgs couplings in , the couplings in eq. (34) “collectively” break an $`SU(3)`$ symmetry protecting the Higgs mass. This is necessary for the cancelation of quadratic divergences to occur, and this is why we obtain a top Yukawa coupling which does not require fine tuning of the Higgs mass. The Yukawa couplings $`\lambda _i`$ can be chosen real by redefining the phases of the $`\chi ^c`$. Expanding to first order in the Higgs $`h`$

$`_{top}`$ $`=`$ $`f(\lambda _1\chi _1^c+\lambda _2\chi _2^c)\chi +{\displaystyle \frac{i}{\sqrt{2}}}(\lambda _1\chi _1^c\lambda _2\chi _2^c)h\left(\begin{array}{c}t\hfill \\ b\hfill \end{array}\right)+\mathrm{}`$ (37)
$`=`$ $`m_\chi \chi ^c\chi i\lambda _tt^ch\left(\begin{array}{c}t\hfill \\ b\hfill \end{array}\right)+\mathrm{},`$ (40)
where in the last step we diagonalized the mass matrix by finding the heavy and light linear combinations of the $`\chi ^c`$. We find a $`\chi `$ mass $`m_\chi =\sqrt{\lambda _1^2+\lambda _2^2}f`$, and a top Yukawa coupling, $`\lambda _t=\sqrt{2}\lambda _1\lambda _2/\sqrt{\lambda _1^2+\lambda _2^2}`$. To obtain a sufficiently large top mass both couplings $`\lambda _i`$ must be of order one.

The absence of quadratic divergences to the Higgs mass at one loop from these interactions is again most easily understood by examining Feynman diagrams with external $`\mathrm{\Phi }`$’s. The diagram in Figure 4.a is quadratically divergent but preserves both $`SU(3)`$ symmetries. Thus it does not contribute to PNGB masses. The diagram 4.b does contribute to the Higgs mass, but it is only log divergent, contributing to the operator $`\lambda _1^2\lambda _2^2/16\pi ^2|\mathrm{\Phi }_2^{}\mathrm{\Phi }_1|^2`$ which contains a Higgs mass of order $`f\lambda /4\pi `$.

Alternatively, one can perform this computation after expanding the $`\mathrm{\Phi }`$’s. In “component form” the vanishing of the quadratic divergence involves “miraculous” cancellations between loops of top quarks and loops of $`\chi `$’s. We demonstrate this calculation for the simplifying choice of $`\lambda _1=\lambda _2\lambda _t`$. Expanded out to the relevant order, the Lagrangian contains

$$\lambda _t(\sqrt{2}f\frac{1}{\sqrt{2}f}h^{}h)\chi ^c\chi +\lambda _tt^ch\left(\begin{array}{c}t\hfill \\ b\hfill \end{array}\right).$$
(41)
In addition to the usual top loop there is also a quadratically divergent $`\chi `$ loop (Figure 5). In the diagram the $`\lambda _t\sqrt{2}f`$ mass insertion on the $`\chi `$ line combines with the $`\lambda _t/\sqrt{2}f`$ from the vertex to exactly cancel the quadratic divergence from the top loop.

Note that although the cancellation involves a higher dimensional operator, its coefficient is related to the Yukawa coupling by the non-linearly realized axial $`SU(3)`$ symmetry.

### 2.2 $`[SU(3)/SU(2)]^2`$ symmetry

Here we give an elegant way of understanding the absence of quadratic divergences in our theory which relies on the structure of the explicit breaking of the $`[SU(3)]^2`$ symmetry acting on the $`\mathrm{\Phi }`$ fields. The $`[SU(3)]^2`$ symmetry not only protects both sets of NGBs ($`\mathrm{\Theta }_{eaten}`$ and $`\mathrm{\Theta }`$) from obtaining a mass, it also forbids any non-derivative couplings of the NGBs. Therefore, in order to generate gauge- and Yukawa couplings for the Higgs, these symmetries must be explicitly broken. The trick is to do this breaking “non-locally in theory space”  or “collectively”. By collective breaking we mean that no single coupling in the Lagrangian breaks the symmetry by itself. Always at least two couplings are required to break any of the two $`SU(3)`$’s.

To see this explicitly, consider the gauge interactions

$$\left|(_\mu +ig_1A_\mu )\mathrm{\Phi }_1\right|^2+\left|(_\mu +ig_2A_\mu )\mathrm{\Phi }_2\right|^2,$$
(42)
where for clarity we have labeled the gauge couplings of $`\mathrm{\Phi }_1`$ and $`\mathrm{\Phi }_2`$ differently. Of course, gauge invariance requires them to have the same value, but as spurions they break different symmetries, and it is useful to keep track of them independently.

The key is to note that if either of the two couplings $`g_1`$ and $`g_2`$ is set to zero the Lagrangian has an exact $`[SU(3)]^2`$ symmetry, and therefore two sets of exact NGBs ($`\mathrm{\Theta }_{eaten}`$ and $`\mathrm{\Theta }`$) result from the spontaneous symmetry breaking. Only with both $`g`$’s non-vanishing are the symmetries explicitly broken to the diagonal gauged $`SU(3)`$. For $`g_2=0`$ we have the two independent symmetries

$`\mathrm{\Phi }_1U_1\mathrm{\Phi }_1,A_\mu U_1A_\mu U_1^{},\mathrm{\Phi }_2U_2\mathrm{\Phi }_2,`$ (43)
while for $`g_1=0`$ we have the symmetries

$`\mathrm{\Phi }_1U_1\mathrm{\Phi }_1,A_\mu U_2A_\mu U_2^{},\mathrm{\Phi }_2U_2\mathrm{\Phi }_2.`$ (44)
Thus when either of the $`g_i`$ is set to zero $`\mathrm{\Theta }`$ is an exact NGB. Any loop correction to its mass must be proportional to both $`g_1`$ and $`g_2`$. But, there are no quadratically divergent one loop diagrams involving both $`g_1`$ and $`g_2`$.

The argument for the Yukawa couplings is very similar. In absence of either $`\lambda _1`$ or $`\lambda _2`$ the Yukawa couplings in eq. (34) preserve two $`SU(3)`$ symmetries. For example, when $`\lambda _2=0`$ we have

$`\mathrm{\Phi }_1U_1\mathrm{\Phi }_1,QU_1Q,\mathrm{\Phi }_2U_2\mathrm{\Phi }_2,`$ (45)
where $`Q`$ denotes the quark triplet. Thus any $`\mathrm{\Theta }`$ mass which is generated from Yukawa loops must be proportional to both $`\lambda _1`$ and $`\lambda _2`$. In fact, it must be proportional to $`|\lambda _1|^2|\lambda _2|^2`$. This follows from the two spurious $`U(1)_i`$ symmetries under which only $`\lambda _i`$ and $`\chi _i^c`$ are charged such that $`\lambda _i\chi _i^c`$ are neutral. Any operator which is generated by loops must also respect these spurious symmetries, and if it doesn’t contain $`\chi `$ s it can only depend on the invariant combinations $`|\lambda _1|^2`$ and $`|\lambda _2|^2`$. Again, there is no quadratically divergent one-loop diagram proportional to $`|\lambda _1|^2|\lambda _2|^2`$.

### 2.3 Lets make the Standard Model

Three ingredients are still missing to turn the above model into a fully realistic extension of the Standard Model. We need i. hypercharge and color, ii. Yukawa couplings for the other fermions, and iii. a quartic self coupling for the Higgs in order to stabilize its VEV.

i. Hypercharge and color: Adding color is trivial. To one loop the only relevance of color for the Higgs mass is a color factor of three in the top quark loop. Adding hypercharge is also straightforward. We introduce it by gauging a $`U(1)_X`$ symmetry which commutes with the $`SU(3)`$ and under which the $`\mathrm{\Phi }_i`$ have charge $`1/3`$. In the non-linear sigma model $`U(1)_X`$ is non-linearly realized but the linear combination corresponding to hypercharge

$`Y={\displaystyle \frac{1}{\sqrt{3}}}T_8+{\displaystyle \frac{1}{3}}X`$ (46)
is unbroken. Here $`T_8=\frac{1}{\sqrt{3}}diag(\frac{1}{2},\frac{1}{2},1)`$ is one of the broken $`SU(3)`$ generators. The $`U(1)_X`$ charges of the fermions are uniquely determined from their hypercharges. For example, $`\chi _i^c`$ has charge $`\frac{2}{3}`$ and $`Q`$ has charge $`\frac{1}{3}`$.

Gauging $`U(1)_X`$ does not introduce new quadratic divergences for the Higgs mass. This is because $`U(1)_X`$ commutes with all the spurious global $`SU(3)`$ symmetries which we used to argue for the absence of divergences. Thus these arguments go through unchanged.

ii. Yukawa couplings for light fermions: Since the Yukawa couplings of the other fermions are small, no care needs to be taken in their coupling to the Higgs in order to avoid large Higgs mass corrections. One possibility is to enlarge all left-handed fermion doublets into triplets and then write Yukawa couplings for up-type quarks and neutrinos just as we did for the top quark. For the bottom quark we can write

$`{\displaystyle \frac{\lambda _b}{\mathrm{\Lambda }}}b^cϵ_{ijk}\mathrm{\Phi }_1^i\mathrm{\Phi }_2^jQ^k+h.c.,`$ (47)
and analogous terms for all down-type quarks and charged leptons. We postpone a more detailed discussion of flavor until Section 4.

iii. Higgs quartic coupling: One possibility for the quartic coupling is to ignore the numerically not very significant one-loop divergence from the Higgs loop and add a quartic coupling without canceling its divergence. A more satisfying model in which this divergence is also canceled requires more work. It is easy to appreciate the difficulty by looking at possible terms which one might add to the Lagrangian to generate the quartic coupling. This would-be Higgs potential is an arbitrary gauge invariant polynomial in $`\mathrm{\Phi }_1`$ and $`\mathrm{\Phi }_2`$. The only non-trivial gauge invariant contraction at our disposal is $`\mathrm{\Phi }_1^{}\mathrm{\Phi }_2`$. The others either vanish ($`ϵ_{ijk}\mathrm{\Phi }_1^i\mathrm{\Phi }_1^j\mathrm{\Phi }_2^k=0`$) or are constant ($`\mathrm{\Phi }_i^{}\mathrm{\Phi }_i=f^2`$). Thus the potential is a function of

$`\mathrm{\Phi }_1^{}\mathrm{\Phi }_2=f^2+if\eta h^{}h{\displaystyle \frac{1}{2}}\eta ^2+\mathrm{}+{\displaystyle \frac{1}{6f^2}}(h^{}h)^2+\mathrm{}`$ (48)
and it’s hermitian conjugate. Focusing on the $`h`$-dependence, we see that the quartic coupling always comes accompanied by a mass term when expanding out a term like $`(\mathrm{\Phi }_1^{}\mathrm{\Phi }_2)^n`$. Setting the coefficient of the quartic to one, the mass is of order $`f`$ which is much too large. Note that if the constant term in the expansion of $`\mathrm{\Phi }_1^{}\mathrm{\Phi }_2`$ were not there, one could generate a quartic by simply squaring eq. (48). This observation will be the key to constructing a quartic in Section 4. Of course, it is possible to fine tune the mass term away. For example the potential $`|\mathrm{\Phi }_1^{}\mathrm{\Phi }_2f^2|^2`$ does not contain a mass term. But this is no better than the tuning of the Higgs mass in the Standard Model as the relative size of coefficients in $`|\mathrm{\Phi }_1^{}\mathrm{\Phi }_2f^2|^2`$ is not stable under quantum corrections.

## 3 A Complete Model: SU(4)

Now we present a complete model of electroweak symmetry breaking with no fine-tuning required to separate the Higgs mass from the cutoff. The main distinction from the previous model is the extension of the gauge group (and therefore the approximate global symmetries) from $`SU(3)`$ to $`SU(4)`$. The expansion of the group allows one to use a new mechanism to generate a quartic coupling.

Take $`\mathrm{\Phi }`$ and $`\mathrm{\Psi }`$ as fields in the non-linear sigma model $`[SU(4)/SU(3)]^2`$ with the diagonal $`SU(4)`$ gauged. The important distinction from the previous model is that the $`SU(4)`$-breaking is not aligned,

$`\mathrm{\Phi }=e^{i\phi /f}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ f\hfill \\ 0\hfill \end{array}\right)`$ (53)
$`\mathrm{\Psi }=e^{i\psi /f}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ 0\hfill \\ f\hfill \end{array}\right)`$ (58)
and only the gauged $`SU(2)`$ is linearly realized. Note that the product $`\mathrm{\Phi }^{}\mathrm{\Psi }`$ contains no constant term! This is exactly the success we were attempting to achieve in the $`SU(3)`$ model. Raised to the appropriate power, it could potentially contain a term quartic in Higgses without a quadratic term.

This (mis-)alignment is stabilized when one adds the interaction $`|\mathrm{\Phi }^{}\mathrm{\Psi }|^2`$ to the potential with a positive coefficient. As we will see, this interaction becomes a positive squared mass for a charged scalar field – the only uneaten (complex) scalar in this example. The two $`SU(2)`$ doublets which live in $`\mathrm{\Phi }`$ and $`\mathrm{\Psi }`$ remain exact NGBs (they are in fact eaten) and thus this term does not induce a quartic term.

To reproduce the successes of the $`SU(3)`$ model (and produce uneaten Higgs doublets), we break the gauged $`SU(4)SU(2)`$ twice. Our model has four sets of sigma model fields, ($`\mathrm{\Phi }_i,\mathrm{\Psi }_i,i=1,2`$). Each contains one complex $`SU(2)`$ doublet. Of the four doublets, two are eaten by the heavy $`SU(4)`$ gauge bosons, leaving a two-Higgs-doublet model.

The complete counting goes as follows: the $`[SU(4)/SU(3)]^4`$ represents $`(158)\times 4=28`$ real components, 12 of which are eaten when the $`SU(4)`$ gauge group is broken to $`SU(2)`$. The remaining 16 consist of two complex doublets $`h_u`$ and $`h_d`$, three complex $`SU(2)`$ singlets $`\sigma _1`$, $`\sigma _2`$ and $`\sigma _3`$, and two real scalars $`\eta _u`$ and $`\eta _d`$. One possible parameterization is as follows:

$`\mathrm{\Phi }_1=e^{+i/f}e^{+i\mathrm{\Sigma }_1/f}e^{+i\mathrm{\Sigma }_2/f}e^{+i\mathrm{\Sigma }_3/f}e^{+i\eta _u/f}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ f\hfill \\ 0\hfill \end{array}\right)`$ (63)
$`\mathrm{\Phi }_2=e^{i/f}e^{+i\mathrm{\Sigma }_1/f}e^{i\mathrm{\Sigma }_2/f}e^{i\mathrm{\Sigma }_3/f}e^{i\eta _u/f}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ f\hfill \\ 0\hfill \end{array}\right)`$ (68)
$`\mathrm{\Psi }_1=e^{+i/f}e^{i\mathrm{\Sigma }_1/f}e^{+i\mathrm{\Sigma }_2/f}e^{i\mathrm{\Sigma }_3/f}e^{+i\eta _d/f}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ 0\hfill \\ f\hfill \end{array}\right)`$ (73)
$`\mathrm{\Psi }_2=e^{i/f}e^{i\mathrm{\Sigma }_1/f}e^{i\mathrm{\Sigma }_2/f}e^{+i\mathrm{\Sigma }_3/f}e^{i\eta _d/f}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ 0\hfill \\ f\hfill \end{array}\right)`$ (78)
with

$`={\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{cccc}\begin{array}{cc}0\hfill & 0\hfill \\ 0\hfill & 0\hfill \end{array}& h_u& h_d& \\ h_u^{}& 0& 0& \\ h_d^{}& 0& 0& \end{array}\right)`$ $`\mathrm{\Sigma }_1={\displaystyle \frac{1}{2}}\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& \sigma _1\\ 0& 0& \sigma _1^{}& 0\end{array}\right)`$ (88)
$`\mathrm{\Sigma }_2={\displaystyle \frac{1}{2}}\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& \sigma _2\\ 0& 0& \sigma _2^{}& 0\end{array}\right)`$ $`\mathrm{\Sigma }_3={\displaystyle \frac{1}{2}}\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& \sigma _3\\ 0& 0& \sigma _3^{}& 0\end{array}\right)`$ (97)
$`\eta _u={\displaystyle \frac{\eta _u}{6}}\left(\begin{array}{cccc}\hfill 1& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 3& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1\end{array}\right)`$ $`\eta _d={\displaystyle \frac{\eta _d}{6}}\left(\begin{array}{cccc}\hfill 1& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 1& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 3\end{array}\right).`$ (106)
A quartic term for this model will come from the four couplings

$`\kappa _{11}|\mathrm{\Phi }_1^{}\mathrm{\Psi }_1|^2+\kappa _{22}|\mathrm{\Phi }_2^{}\mathrm{\Psi }_2|^2+\kappa _{12}|\mathrm{\Phi }_1^{}\mathrm{\Psi }_2|^2+\kappa _{21}|\mathrm{\Phi }_2^{}\mathrm{\Psi }_1|^2.`$ (107)
To leading order in each field, these operators produce the potential:

$`\kappa _{11}f^2|\sigma _1+\sigma _3|^2`$ $`+`$ $`\kappa _{22}f^2|\sigma _1\sigma _3|^2`$
$`+\kappa _{12}f^2|\sigma _1+\sigma _2ih_u^{}h_d/f|^2`$ $`+`$ $`\kappa _{21}f^2|\sigma _1\sigma _2ih_u^{}h_d/f|^2.`$
Any one coupling does not produce a potential for the Higgses. In fact, removing any of the couplings above removes the quartic Higgs term. To see this, note that there are other parameterizations in which the $`h_u^{}h_d`$ appear in different operators.

These couplings generate masses of order $`\sqrt{\kappa }f`$ for the complex scalars $`\sigma _i`$ as well as trilinear terms marrying these scalars to $`h_u^{}h_d`$. Integrating out the singlets produces a quartic coupling:

$$\lambda =\frac{4}{\kappa _{11}^1+\kappa _{22}^1+\kappa _{12}^1+\kappa _{21}^1}$$
(108)
which vanishes as any one $`\kappa _{ij}0`$, thus all four terms are required to produce a tree-level quartic term.

From symmetry arguments we see why this works. The non-linearly realized global symmetry of the model is approximately $`[SU(4)/SU(3)]^4`$ and contains, among other things, four Higgs-like doublets which are NGBs. Two of the doublets are eaten. A single operator, e.g., $`|\mathrm{\Phi }_1^{}\mathrm{\Psi }_1|^2`$, explicitly breaks the symmetry down to $`[SU(4)/SU(3)]^2\times [SU(4)/SU(2)]`$, but this symmetry also contains four doublet NGBs. Adding $`|\mathrm{\Phi }_2^{}\mathrm{\Psi }_2|^2`$ leaves $`[SU(4)/SU(2)]^2`$ again producing four doublets. The existence of three of the four operators breaks enough symmetry to allow a quartic term, but only a small one is induced at loop level.

This structure also suppresses one-loop contributions to the Higgs mass. The symmetry arguments above show that more than one spurion is required to generate a Higgs potential, and therefore a Higgs mass. Thus, at one loop there are no quadratically divergent diagrams.

The gauge and top loops work similarly to those in the $`SU(3)`$ model. The gauge loop is canceled for each Higgs doublet by the massive gauge bosons and the calculation goes through as in Section 2. Gauging an additional $`U(1)`$ symmetry results in the existence of hypercharge at the weak scale and does not contribute to the Higgs potential. The top Yukawa coupling could, for example, come from:

$`_{top}=(\lambda _1\chi _1^c\mathrm{\Phi }_1^{}+\lambda _2\chi _2^c\mathrm{\Phi }_2^{}+\lambda _3\chi _3^c\mathrm{\Psi }_1^{})Q`$ (109)
where $`Q^T=(t,b,\chi _1,\chi _2)`$. The $`\chi _3^c`$ field is there to cancel the hypercharge anomaly and to marry and give a mass to $`\chi _2`$. The $`\lambda _3`$ term does not play a role in generating a top Yukawa coupling and thus the coupling can be taken to the cutoff (i.e., $`\lambda _34\pi `$) thus decoupling these extra fields and making the physics the same as the $`SU(3)`$ case of the previous section. The bottom Yukawa coupling can be written as

$`_{bottom}=\lambda _bb^cϵ_{ijkl}\mathrm{\Phi }_1^i\mathrm{\Psi }_1^j\mathrm{\Psi }_2^kQ^l.`$ (110)
This single coupling generates a quadratic divergence to the down-type Higgs mass of the form $`|\mathrm{\Psi }_1^{}\mathrm{\Psi }_2|^2`$. For a small bottom Yukawa coupling ($`\lambda _b\stackrel{<}{_{}}0.1`$), the contribution remains at or below the weak scale.

### 3.1 Electroweak Symmetry Breaking

The purpose of the quartic term in the Higgs potential is to stabilize the Higgs VEV . The $`SU(4)`$ model above has a quartic potential of the form

$$_{quartic}=\lambda |h_u^{}h_d|^2$$
(111)
where $`h_u`$ and $`h_d`$ are $`SU(2)`$ doublets with hypercharge $`Y=1/2`$. Successful electroweak symmetry breaking requires the existence of a mass term of the form $`Bh_u^{}h_d`$ where $`B`$ is of order $`M_W^2`$. Such a term comes from operators $`B_{11}\mathrm{\Phi }_1^{}\mathrm{\Psi }_1`$, $`B_{22}\mathrm{\Phi }_2^{}\mathrm{\Psi }_2`$, $`B_{12}\mathrm{\Phi }_1^{}\mathrm{\Psi }_2`$, and $`B_{21}\mathrm{\Phi }_2^{}\mathrm{\Psi }_1`$, and can be written with the quartics as $`\kappa _{11}|b_{11}+\mathrm{\Phi }_1^{}\mathrm{\Psi }_1|^2+\kappa _{22}|b_{22}+\mathrm{\Phi }_2^{}\mathrm{\Psi }_2|^2+\mathrm{}`$. These new operators contain linear terms for the uncharged scalar fields $`\sigma _i`$ causing them to obtain VEVs. When the $`\sigma `$’s are shifted to their minima, a mass term of the form $`Bh_u^{}h_d`$ is produced with

$$B=\frac{1}{2}\lambda \underset{ij}{}b_{ij}.$$
(112)
Here – for simplicity – we have taken the $`b_i`$ to be real, though in general their phases may have interesting implications for CP violation. These $`B`$ terms are themselves spurions which explicitly break more symmetries than the quartic couplings. To see this, note that they are the only terms so far which produce a potential for the $`\eta `$ fields. Thus their size, which needs to be of order $`f^2/16\pi ^2`$ is technically natural, but undetermined from dynamics in the effective theory below $`\mathrm{\Lambda }`$.

In addition, operators of the form $`|\mathrm{\Phi }_1^{}\mathrm{\Phi }_2|^2`$ and $`|\mathrm{\Psi }_1^{}\mathrm{\Psi }_2|^2`$ are generated by two-loop quadratic-divergent and one-loop log-divergent diagrams. They produce the mass terms

$$_{mass}=m_2^2|h_u|^2+m_1^2|h_d|^2$$
(113)
with a natural size of order $`f^2/16\pi ^2`$. We require $`m_2^2,m_1^2>0`$ (or else the Higgs vev could run away to $`f`$). Electroweak symmetry breaking occurs if

$$B>\sqrt{m_2^2m_1^2}.$$
(114)
This Higgs potential is of the same form as the one in the $`SU(6)/Sp(6)`$ little Higgs model , and we repeat some of the phenomenology here. Formulas for the scalar masses in general two Higgs doublet models are conveniently collected in . Minimizing the potential under these conditions gives $`\mathrm{tan}\beta =v_u/v_dh_u/h_d=\sqrt{m_1^2/m_2^2}`$ and

$$\frac{2B}{\mathrm{sin}2\beta }=m_2^2+m_1^2+2\lambda v^2$$
(115)
where $`v=174`$ GeV is the electroweak symmetry breaking scale.

The masses of the two CP-even Higgs bosons are

$$M_{h^0}^2,M_{H^0}^2=\frac{B}{\mathrm{sin}2\beta }\pm \sqrt{\frac{B^2}{\mathrm{sin}^22\beta }+\lambda v^2\left(\lambda v^2\frac{2B}{\mathrm{sin}2\beta }\right)\mathrm{sin}^22\beta }.$$
(116)
The lightest CP-even Higgs boson is bounded from above by $`M_{h^0}^2\lambda v^2`$. This bound is saturated for $`m_1^2=m_2^2\mathrm{sin}2\beta =1`$. The CP-odd and charged Higgses have masses

$`M_{A^0}^2`$ $`=`$ $`{\displaystyle \frac{2B}{\mathrm{sin}2\beta }}`$ (117)
$`M_{H^\pm }^2`$ $`=`$ $`M_{A^0}^2\lambda v^2`$ (118)
### 3.2 The Standard Model Embedding

Now we can construct a complete standard model based on the $`SU(4)`$ theory. Collecting the pieces together, the Lagrangian of the theory is

$$=_{kinetic}+_{quarks}+_{leptons}+_{Higgs}$$
(119)
with

$$_{kinetic}=|𝒟_\mu \mathrm{\Phi }_i|^2+|𝒟_\mu \mathrm{\Psi }_i|^2+\left[\mathrm{fermion}\mathrm{and}\mathrm{gauge}\mathrm{kinetic}\mathrm{terms}\right].$$
(120)
Here $`𝒟_\mu =(_\mu +ig_4A_\mu ^aT^ai\frac{g_X}{4}A_\mu ^X)`$ and the $`\frac{1}{4}`$ in the coupling of $`A_\mu ^X`$ represents the $`U(1)_X`$ charge of $`\mathrm{\Phi }`$ and $`\mathrm{\Psi }`$. $`A_\mu `$ and $`A_\mu ^X`$ are the $`SU(4)`$ and $`U(1)_X`$ gauge fields respectively. The Yukawa couplings for quarks appear as

$$_{quarks}=(\lambda _1^u\chi _{u1}^c\mathrm{\Phi }_1^{}+\lambda _2^u\chi _{u2}^c\mathrm{\Phi }_2^{}+\lambda _3^u\chi _{u3}^c\mathrm{\Psi }_1^{})Q+\lambda ^dd^c\mathrm{\Phi }_1\mathrm{\Psi }_1\mathrm{\Psi }_2Q$$
(121)
with $`Q=(q,\chi _{u1},\chi _{u2})^T`$. We have suppressed flavor and $`SU(4)`$ indices for clarity. The $`\lambda `$ couplings are $`3\times 3`$ matrices in flavor space – the combination of the first three produces the standard Yukawa matrix while $`\lambda ^d`$ is simply the Yukawa matrix for down-type quarks. Similarly, for the charged leptons:

$$_{leptons}=(\lambda _1^\nu \chi _{\nu 1}^c\mathrm{\Phi }_1^{}+\lambda _2^\nu \chi _{\nu 2}^c\mathrm{\Psi }_1^{})Q+\lambda ^ee^c\mathrm{\Phi }_1\mathrm{\Psi }_1\mathrm{\Psi }_2L$$
(122)
where $`L=(\mathrm{},\chi _{\nu 1},\chi _{\nu 2})^T`$, we will discuss neutrino masses in the next section. Finally, the tree level scalar potential is

$`_{scalar}={\displaystyle \underset{ij}{}}\kappa _{ij}|b_{ij}+\mathrm{\Phi }_i^{}\mathrm{\Psi }_j|^2`$ (123)
as discussed in the previous subsection.

Hypercharge is a linear combination of the $`SU(4)`$ generator

$$T^{15}=\sqrt{2}diag(1/4,1/4,+1/4,+1/4)$$
(124)
and the external $`U(1)_X`$. Thus the $`X`$ charges of the $`SU(4)`$-singlet fermions are just their respective hypercharges while the $`X`$ charges of $`SU(4)`$ vectors are the hypercharges of the $`SU(2)`$ doublets they contain plus $`1/4`$. Explicitly, $`(\mathrm{\Phi }_i,\mathrm{\Psi }_i,L,Q)`$ have $`X`$ charges $`(1/4,1/4,1/4,+5/12)`$.

## 4 Spectrum and Constraints

A complete analysis of the phenomenology is beyond the scope of this paper but we would like to report on our initial explorations in this direction. The results are encouraging: we find significant constraints but there are large regions of parameter space which are in agreement with experiment while at the same time solving the hierarchy problem. Interestingly, the preferred region or parameter space will be directly explored at the Tevatron and LHC. More specifically, we will discuss i. precision electroweak constraints ii. the spectrum and direct searches iii. flavor physics.

i. precision constraints: One of the most stringent constraints on models of new physics at the TeV scale comes from isospin violating couplings of the light fermions to the $`W`$ and $`Z`$ bosons. One source of isospin violation is the different treatment of the Yukawa couplings for up and down-type quarks in eq. (121). Through the Higgs vev up-type quarks mix with the heavy $`\chi `$ fermions whereas down-type quarks don’t. As we will now show, this mixing prefers different scales $`f_i`$ for the different $`\mathrm{\Phi }_i`$.

For simplicity, we revert to our $`SU(3)`$ model where the same mixing occurs. The up-type Yukawa couplings are

$$(\lambda _1u_1^c\mathrm{\Phi }_1^{}+\lambda _2u_2^c\mathrm{\Phi }_2^{})\left(\begin{array}{c}u\hfill \\ d\hfill \\ \chi \hfill \end{array}\right),$$
(125)
where $`\lambda _i`$ are $`3\times 3`$ matrices in flavor space. In order to avoid large flavor changing effects (see FCNC discussion below) we take $`\lambda _2`$ proportional to the unit matrix and of order one whereas $`\lambda _1`$ is approximately equal to the usual Yukawa couplings of the up-type quarks in the Standard Model. Furthermore, we allow different scales $`f_1`$ and $`f_2`$ for the non-linear sigma models. The heavy $`SU(3)`$ gauge bosons eat a linear combination of the NGBs which resides mostly in the sigma model with the larger scale, and the little Higgs lives mostly in the sigma model with the smaller scale:

$`\mathrm{\Phi }_1=e^{i\mathrm{\Theta }\frac{f_2}{f_1}}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ f_1\hfill \end{array}\right),\mathrm{\Phi }_2=e^{i\mathrm{\Theta }\frac{f_1}{f_2}}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ f_2\hfill \end{array}\right)`$ (132)
where

$`\mathrm{\Theta }=\left(\begin{array}{cc}\begin{array}{cc}0\hfill & 0\hfill \\ 0\hfill & 0\hfill \end{array}& h\\ h^{}& 0\end{array}\right)/f_{12}\mathrm{and}f_{12}^2=f_1^2+f_2^2.`$ (137)
Substituting the Higgs by its expectation value $`h^T=(v,0)`$ we obtain the mass matrix

$`\left(u_1^cu_2^c\right)\left(\begin{array}{cc}\lambda _1v\frac{f_2}{f_{12}}& \lambda _1f_1\\ \lambda _2v\frac{f_1}{f_{12}}& \lambda _2f_2\end{array}\right)\left(\begin{array}{c}u\hfill \\ \chi \hfill \end{array}\right)`$ (142)
Since $`\lambda _2>>\lambda _1`$ for the light quarks we see that the heavy quarks are approximately $`u_2^c,\chi `$ with masses $`\lambda _2f_21`$ TeV, and the light (SM) quarks are $`u_1^c,u`$ with masses $`\lambda _1vf_2/f_{12}`$. In addition, there is small mixing between light and heavy quarks. The mixing between the $`u^c`$ fields is not physical and can be removed by a change of basis. However, mixing between the $`SU(2)`$ doublet component $`u`$ and the singlet $`\chi `$ is significant because it alters the couplings of up-type quarks to the $`W`$ and $`Z`$. The mixing angle is $`vf_1/(f_2f_{12})`$, which reduces the coupling of an up-type quark by

$`\delta g={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{f_1v}{f_2f_{12}}}\right)^2.`$ (143)
For $`f_1f_21`$ TeV and $`v=175`$ GeV the shift in the coupling is 1%. A similar shift also occurs in the couplings of neutrinos from their mixing with heavy partners. This is problematic because precision measurements at LEP and SLC have determined the gauge couplings of light fermions to a precision of $`2\times 10^3`$ . However, we also see that it is easy to strongly suppress the mixing by taking $`f_2>f_1`$. For example, taking $`f_2=2`$ TeV and $`f_1=1`$ TeV we have $`\delta g10^3`$. We see that the part of parameter space with $`f_2>f_1`$ is preferred.

We should check that taking unequal $`f_i`$ allows a large enough top Yukawa coupling and does not destabilize the Higgs mass. Diagonalizing eq. (142) to leading order in $`v^2/f^2`$ for the third generation we obtain the mass of the heavy partner of the top and the top Yukawa coupling

$`m_\chi =\sqrt{\lambda _1^2f_1^2+\lambda _2^2f_2^2},\lambda _t=\lambda _1\lambda _2\sqrt{{\displaystyle \frac{f_1^2+f_2^2}{\lambda _1^2f_1^2+\lambda _2^2f_2^2}}}`$ (144)
Happily, a wide range of $`f_i`$ and $`\lambda _i`$ can give a heavy top quark without fine tuning of the Higgs mass. For example, taking $`f_1=.5`$ TeV, $`f_2=2`$ TeV, $`\lambda _1=\sqrt{2}`$ and $`\lambda _2=1/3`$ we obtain the correct top Yukawa. To estimate the degree of fine-tuning recall that the top loop contribution to the Higgs mass is cut off by $`m_\chi `$. Thus $`\delta m_h^2m_\chi ^2\lambda _t^2/16\pi ^2`$ which requires no fine tuning for $`m_\chi 1`$ TeV.

We now turn to computing the masses and mixings of the gauge bosons in the full $`SU(4)\times U(1)_X`$ model. Transitions mediated by the heavy $`SU(4)`$ gauge bosons contribute to precision electroweak measurements leading to constraints on the $`f_i`$. A useful parametrization of the non-linear sigma model fields $`\mathrm{\Phi }_i`$ and $`\mathrm{\Psi }_i`$ with general $`f_i`$ is

$`\mathrm{\Phi }_1=e^{+i_u\frac{f_2}{f_1}}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ f_1\hfill \\ 0\hfill \end{array}\right)\mathrm{\Phi }_2=e^{i_u\frac{f_1}{f_2}}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ f_2\hfill \\ 0\hfill \end{array}\right)`$ (153)
$`\mathrm{\Psi }_1=e^{+i_d\frac{f_4}{f_3}}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ 0\hfill \\ f_3\hfill \end{array}\right)\mathrm{\Psi }_2=e^{i_d\frac{f_3}{f_4}}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ 0\hfill \\ f_4\hfill \end{array}\right)`$ (162)
where

$`_u=\left(\begin{array}{ccc}\begin{array}{cc}0\hfill & 0\hfill \\ 0\hfill & 0\hfill \end{array}& h_u& \begin{array}{c}0\hfill \\ 0\hfill \end{array}\hfill \\ h_u^{}& 0& 0\hfill \\ \begin{array}{cc}0& 0\end{array}& 0& 0\hfill \end{array}\right)/f_{12}`$ $`_d=\left(\begin{array}{ccc}\begin{array}{cc}0\hfill & 0\hfill \\ 0\hfill & 0\hfill \end{array}& \begin{array}{c}0\hfill \\ 0\hfill \end{array}& h_d\hfill \\ \begin{array}{cc}0& 0\end{array}& 0& 0\hfill \\ h_d^{}& 0& 0\hfill \end{array}\right)/f_{34}`$ (179)
Here we ignore the small contributions to masses from small vevs for the fields $`\sigma _i`$. The photon and the $`Z`$ are linear combinations of four neutral gauge bosons: three gauge bosons which correspond to the diagonal $`SU(4)`$ generators $`T^3=\frac{1}{2}diag(1,1,0,0)`$, $`T^{12}=\frac{1}{2}diag(0,0,1,1)`$, and $`T^{15}=\frac{1}{\sqrt{8}}diag(1,1,1,1)`$, and the $`U(1)_X`$ gauge field $`B_\mu ^x`$. Two linear combinations obtain masses of order $`f`$ from the kinetic terms of the $`\mathrm{\Phi }_i`$ and $`\mathrm{\Psi }_i`$

$`\left|{\displaystyle \frac{g}{2}}A_\mu ^{12}+{\displaystyle \frac{g}{2\sqrt{2}}}A_\mu ^{15}{\displaystyle \frac{g_x}{4}}B_\mu ^x\right|^2f_{12}^2`$
$`\left|{\displaystyle \frac{g}{2}}A_\mu ^{12}+{\displaystyle \frac{g}{2\sqrt{2}}}A_\mu ^{15}{\displaystyle \frac{g_x}{4}}B_\mu ^x\right|^2f_{34}^2`$ (180)
where $`g,g_x`$ are the $`SU(4)`$ and $`U(1)`$ gauge couplings and $`f_{ij}^2=f_i^2+f_j^2`$. In the following we specialize to $`ff_{12}=f_{34}`$ for which the heavy gauge boson mass matrix simplifies. Then the two heavy eigenstates are

$`Z_\mu ^{\prime \prime }`$ $`=`$ $`A_\mu ^{12}`$
$`Z_\mu ^{}`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}gA_\mu ^{15}g_xB_\mu ^x}{\sqrt{2g^2+g_x^2}}},`$
with masses $`m_{Z^{\prime \prime }}=gf`$ and $`m_Z^{}=\frac{gf}{2}\sqrt{2+g_x^2/g^2}`$. The two eigenstates which remain massless at this order are

$`W_\mu ^3`$ $`=`$ $`A_\mu ^3`$
$`B_\mu `$ $`=`$ $`{\displaystyle \frac{g_xA_\mu ^{15}+\sqrt{2}gB_\mu ^x}{\sqrt{2g^2+g_x^2}}}.`$ (181)
The $`Z`$ obtains its mass from the Higgs vevs $`v=\sqrt{v_u^2+v_d^2}`$. Ignoring mixing with the $`Z^{}`$ the mass term is

$`v^2\left|{\displaystyle \frac{g}{2}}W_\mu ^3{\displaystyle \frac{g_x}{2}}{\displaystyle \frac{1}{\sqrt{1+g_x^2/2g^2}}}B_\mu \right|^2.`$ (182)
From this expression we can read off the standard model gauge couplings. We see that the $`SU(2)`$ coupling of the standard model is equal to the $`SU(4)`$ coupling $`g`$ and – setting the coefficient of $`B_\mu `$ equal to $`g^{}/2`$ – we have

$$g^{}=g_x/\sqrt{1+\frac{g_x^2}{2g^2}}.$$
(183)
Deviations from the standard model arise in this model at order $`v^2/f^2`$ from mixing of the $`Z`$ with the $`Z^{}`$. Explicitly, the mixing is determined by diagonalizing the $`Z`$$`Z^{}`$ mass matrix

$$\frac{g^2}{2}(\begin{array}{cc}v^2(1+t^2)& v^2(1t^2)\sqrt{1+t^2}/\sqrt{2t^2}\\ v^2(1t^2)\sqrt{1+t^2}/\sqrt{2t^2}& 2f^2/(2t^2)\end{array})$$
(184)
where $`t=g^{}/g=\mathrm{tan}\theta _W`$ and $`\theta _W`$ is the weak mixing angle. Diagonalizing, we find a contribution to $`\delta \rho `$ from the shift in the $`Z`$ mass

$$\delta \rho \frac{\delta m_W^2}{m_W^2}\frac{\delta m_Z^2}{m_Z^2}=\frac{v^2}{2f^2}(1t^2)^2+1.510^3\left(\frac{2.2\mathrm{TeV}}{f}\right)^2.$$
(185)
Given that a standard model fit predicts a $`W`$-mass which is lower than the experimental value by about $`1.6\sigma `$ , this correction actually improves the precision electroweak fit for $`f2.2`$ TeV. Alternatively, demanding a fit that is at least as good as the standard model implies a bound of $`f\stackrel{>}{_{}}1.5`$ TeV.

Another observable affected by the new gauge bosons are four-fermion operators. The bound on new contributions to the four-electron operator, for example, is quite severe. The exchange of the $`Z^{}`$ produces an operator of the size:

$$\frac{(1t^2)^2}{8f^2}\overline{e}\gamma ^\mu e\overline{e}\gamma _\mu e$$
(186)
for left-left currents. Using current bounds on this operator we find the requirement that $`f\stackrel{>}{_{}}1.5`$ TeV.

ii. the spectrum and direct searches: In the UV the standard model $`SU(2)\times U(1)`$ gauge group is enlarged to $`SU(4)\times U(1)`$. Thus there are 12 new massive gauge bosons with masses near a TeV. Two of them are the $`Z^{}`$ and $`Z^{\prime \prime }`$ discussed above. For the parameter choice $`f_1^2+f_2^2=f_3^2+f_4^2f^2`$ the $`Z^{\prime \prime }`$ has mass $`gf`$ and does not couple to standard model fermions. The $`Z^{}`$ has mass $`.77gf`$ and couples to quarks and leptons. At the Tevatron it would appear as an s-channel resonance which decays to pairs of leptons. The limit on the mass of such a $`Z^{}`$ from CDF  is in the 700-800 GeV range, implying a bound $`f\stackrel{>}{_{}}`$ 1 TeV. The off-diagonal $`SU(4)`$ gauge bosons and their masses are

$`\left(\begin{array}{cccc}& & Y^0& Y^0\\ & & X^{}& X^{}\\ \overline{Y}^0& X^+& & Y^{\prime \prime 0}\\ \overline{Y}^0& X^+& \overline{Y}^{\prime \prime 0}& \end{array}\right)=\left(\begin{array}{cccc}& & .5& .5\\ & & .5& .5\\ .5& .5& & 1\\ .5& .5& 1& \end{array}\right)gf`$ (195)
The $`Y^{\prime \prime 0}`$ only couples to the heavy fermions and is therefore extremely difficult to detect. All others couple to one light and one heavy fermion. They can be produced in association with a heavy fermion or else appear in t-channel diagrams.

There are two vector-like heavy quarks of charge $`2/3`$ for each generation. As we discussed above, one of them mixes with up-type quarks and can be produced singly in t-channel $`W`$ exchange. The LHC reach in this channel can be as large as several TeV . The masses of these quarks are not completely determined because they depend on unknown Yukawa couplings. But flavor constraints suggest that their masses are generation independent and since naturalness requires a partner for the top quark below $``$ 2 TeV we expect at least one set of these quarks to be visible at the LHC. In terms of model parameters, the new quark masses are $`\sqrt{\lambda _1^2f_1^2+\lambda _2^2f_2^2}`$ and $`\lambda _3f_3`$.

In addition, there are also vector-like heavy leptons which mix with the neutrinos. These fermions are impossible to discover directly but their existence might be inferred by missing energy signals or through their mixing with the light neutrinos in precision data. Note also that this mixing reduces neutral currents of neutrinos more than charged currents which might help to explain the NuTeV anomaly  and the slightly reduced invisible width of the $`Z`$.

The scalar spectrum consists of the two Higgs doublets with masses near the weak scale, three complex neutral fields $`\sigma _i`$ with masses $`\kappa f`$, and two real scalar fields $`\eta _i`$ with masses of order the weak scale. The $`\eta _i`$ only couple to heavy fields and in the $`b`$-terms of the Higgs potential, thus they are very difficult to detect despite their relatively low masses.

iii. flavor: The $`SU(4)`$ model has more Yukawa matrices than the standard model. In general these additional sources of flavor violation lead to flavor changing neutral currents. The easiest way to suppress the flavor violation is to assume that the new Yukawa matrices are proportional to the unit matrix. This assumption imposes constraints on the UV completion of the theory, but it is “technically natural” in the effective theory (loops in the effective theory only generate small corrections). Note that this is similar to the assumption of universal soft masses in supersymmetry.

For example, the quark Yukawa couplings are then

$$_{quarks}=(\lambda ^uu_1^c\mathrm{\Phi }_1^{}+I_1u_2^c\mathrm{\Phi }_2^{}+I_2u_3^c\mathrm{\Psi }_1^{})Q+\lambda ^dd^c\mathrm{\Phi }_1\mathrm{\Psi }_1\mathrm{\Psi }_2Q$$
(196)
where $`\lambda ^u`$ and $`\lambda ^d`$ are similar to the usual standard model Yukawa couplings and $`I_1`$ and $`I_2`$ are approximately proportional to the unit matrix (in flavor space). To see that these couplings do not contain dangerous flavor violation, we go to a new basis in which $`\lambda ^u`$ is diagonal. It is convenient to rotate all four components of $`Q`$ in the same way. $`I_1`$ and $`I_2`$ remain unit matrices if $`u_2^c`$ and $`u_3^c`$ are rotated appropriately. Finally, we also diagonalize $`\lambda ^d`$ with a bi-unitary transformation, but this time we transform only the down-type quarks in $`Q`$. In the new basis all Yukawa couplings are diagonal; flavor violation resides only in the gauge couplings of $`W`$, $`X`$, $`Y`$ to quarks and in couplings of multiple Higgses to quarks which arise from expanding out the down Yukawa operator. The latter couplings are small and only appear in loops. The former are also easily shown to be harmless. There are the usual $`W`$ couplings proportional to the CKM matrix in addition to new couplings between one down-type quark, one heavy vector-like up-type quark and the heavy gauge bosons $`X`$ and $`Y`$. These couplings allow box diagrams and penguins which are similar to corresponding standard model diagrams with $`W`$s replaced by $`X`$s or $`Y`$s. The resulting flavor changing neutral currents are suppressed relative to the standard model ones by the large masses of the heavy gauge bosons $`m_W^2/m_X^2`$ and can be ignored.

A similar analysis of the lepton sector shows that there is also no dangerous lepton flavor violation in the low energy theory. Of course, the theory may also contain direct flavor violating four Fermi operators suppressed by the cut-off $`\mathrm{\Lambda }`$. Such operators are constrained by $`K\overline{K}`$ mixing and CP violation. The experimental bounds on such operators therefore imply constraints on the unknown UV-theory above 10 TeV.

Small neutrino masses can be obtained by including a higher dimensional lepton number violating operator $`(\mathrm{\Phi }^{}L)^2`$ with a small coefficient in the effective theory. This operator might arise from a generalization of the see-saw mechanism in the UV completion: supermassive right handed neutrinos coupled to the operator which interpolates $`(\mathrm{\Phi }^{}L)`$ in the UV theory.

## 5 Discussion

We have seen that, in a little Higgs model, embedding $`SU(2)_{weak}`$ into a simple group (such as $`SU(3)`$) is enough to cancel one-loop quadratic divergences from gauge and (perturbatively coupled) fermion loops. The $`SU(3)`$ model only lacks a quartic. One possibility is to simply ignore the relatively insignificant fine-tuning from the Higgs couplings and add the quartic by hand in “component” fields (i.e., not the full $`\mathrm{\Phi }`$). This coupling need not be very large as an additional contribution to the quartic would come from the log-divergent and finite contributions to the effective potential from the top sector below the scale $`f`$. At worst, this introduces of order 10% fine-tuning.

There are at least three different possibilities for ultra-violet completions to our models. One is a linear sigma model with supersymmetry protecting the scalar masses above the multi-TeV scale. The $`SU(3)`$ model would work well in this case as the quartic would be provided by the $`D`$-term and the Higgs would remain massless at tree-level as long as there are more than two triplets. In addition, the group is simple enough to embed in a unifying theory: $`SU(3+n)_{color}\times SU(2+n)_{weak}\times U(1)`$ gives coupling constant predictions extremely close to those of the MSSM, when charges are normalized to embed into $`SU(5+2n)`$, and matter is in complete representations except for a split fundamental and anti-fundamental .

If one wishes to complete the theory above $`4\pi f`$ with a strongly coupled theory, the coset space we’ve used would require something different than a QCD-like model. For example, a gauged $`SU(7)`$ with four fundamentals, one anti-fundamental and one anti-symmetric tensor produces the symmetry-breaking pattern $`SU(4)SU(3)`$ assuming the fundamentals condense with the anti-fundamental. Fermion and quartic interactions would require additional dynamics as in extended-technicolor .

A more interesting possibility would be a linear sigma model completing into another little Higgs theory at a higher scale, $`F10TeV`$. This may be possible in the $`SU(3)`$ theory if the “fundamental” quartic added need not be too large. One example would be $`[SU(7)/SO(7)]^2`$ with $`SU(3)`$ subgroups gauged similar to the “littlest Higgs” . This could produce two light $`SU(3)`$ triplets with fermion couplings causing vacuum misalignment and $`SU(3)`$-breaking at a scale $`fF/4\pi `$. If a model of this type works, it would provide a weakly coupled theory of electroweak symmetry breaking valid up to $`>`$100 TeV.

## 6 Acknowledgements

We thank Nima Arkani-Hamed, Andy Cohen and Sekhar Chivukula for numerous helpful discussions. DK thanks the visitor program of the particle theory group at the Boston University Physics Department for their hospitality, and we thank the Aspen Center for Physics for providing a stimulating environment during the early phase of this work. MS is supported by the DOE grant DE-FG02-90ER-40560 and the Outstanding Junior Investigator Award DE-FG02-91ER40676.


format: text

text:
# Efficient Algorithm On a Non-staggered Mesh for Simulating Rayleigh-Bénard Convection in a Box

## I Introduction

Experiments over the last three decades have discovered many fascinating and poorly understood examples of pattern formation in large-aspect-ratio Rayleigh-Bénard convection . Because of the prominent role that these experiments play in understanding sustained nonequilibrium systems  and because many of the observed phenomena such as spatiotemporal chaos are difficult to analyze mathematically , there is a need to develop computer codes that can simulate these experiments quantitatively so that theory and experiment can be compared with one another. Once validated, such codes can further be used to explore regimes not easily attained by experiment such as low Prandtl number, and to calculate quantities that are difficult to deduce from experimental data, such as mean flows  and fractal dimensions .

The regime of large aspect ratio $`\mathrm{\Gamma }`$ (ratio of horizontal fluid width to fluid depth) poses significant computational challenges. Many numerical degrees of freedom (basis functions or mesh points) are needed to represent the spatial features of the fluid and often the dynamics needs to be studied over long times (many multiples of the horizontal thermal diffusion time $`\tau _h=\mathrm{\Gamma }^2\tau _v`$, where $`\tau _v=d^2/\kappa `$ is the vertical thermal diffusion time defined in terms of the fluid depth $`d`$ and fluid thermal diffusivity $`\kappa `$) to insure that nontransient behavior is being observed. Since the largest time step allowed by numerical stability for explicit or semi-implicit algorithms (the ones most commonly used in Navier-Stokes calculations) is typically $`0.05\tau _v`$ or smaller, simulations in a representative $`\mathrm{\Gamma }=50`$ cell may require $`10^5`$ or more time steps to eliminate a transient and then study the statistically stationary properties of the asymptotic dynamics. The many degrees of freedom and long integration times together with the need to repeat runs for different parameter values imply that efficient algorithms are essential for studying the large-aspect-ratio regime.

Because of these computational challenges, there have been few simulations of three-dimensional Rayleigh-Bénard convection with aspect ratios exceeding 10. Recent calculations with $`\mathrm{\Gamma }`$ as large as 64 have been carried out by Pesch and collaborators, who used a pseudospectral code on serial and parallel computers to study spiral defect chaos, rotating convection, and other problems . However, their code uses periodic boundaries and so can not take into account quantitatively the influence of lateral walls on the bulk dynamics. Arter and Newell  and Tomita and Abe  have carried out simulations in large boxes with thermally insulated no-slip sidewalls, the former in a $`16\times 11.5`$ aspect ratio box, the latter in a $`\mathrm{\Gamma }=18.84`$ square box. Xi et al  have studied the transition to spatiotemporal chaos of a convecting fluid in a $`\mathrm{\Gamma }=60`$ square cell, but with free-slip horizontal boundaries that are difficult to achieve experimentally. Finally, a Caltech-Duke collaboration has recently reported results  obtained with a parallel spectral element code  for aspect ratios up to 30. Their code can treat quantitatively most geometries and lateral boundaries used by experimentalists, including ramps , spoiler fins , and lateral walls of finite thickness and finite thermal conductivity. However, the generality of the spectral element algorithm makes it substantially more expensive to run than algorithms optimized for a simple geometry such as a box or cylinder.

In this paper, we introduce and analyze an efficient semi-implicit finite-difference algorithm for studying incompressible Rayleigh-Bénard convection in a box, with lateral walls that are periodic, thermally insulated, or thermally conducting. The code complements the more flexible spectral element approach  by being more than an order of magnitude more efficient on a serial processor, for a box with these boundary conditions. It is well suited for studying long-time dynamics of small- to moderate-aspect-ratio boxes  ($`\mathrm{\Gamma }20`$), with lateral boundaries that are close to those of many experiments, although not fully quantitatively accurate since the finite thickness and finite thermal diffusivity of the lateral walls is not taken into account.

The main advantages of our algorithm are the simplicity of implementation and its efficiency on a single processor. Its simplicity arises from the use of a single mesh for all field values. (This is called a “non-staggered” or “collocated” mesh in contrast to a “staggered” mesh for which the values of different fields appear at different points in space ). A non-staggered mesh reduces the effort to write and to validate a code (compared to one using a staggered mesh), and facilitates porting the code to a distributed-memory parallel computer. Earlier work on Navier-Stokes integrators has suggested that non-staggered mesh codes can be numerically unstable because of pressure oscillations . Our results below show that an algorithm to integrate the Boussinesq equations on a non-staggered mesh can be numerically stable.

The use of a single non-staggered mesh also helps to explain the efficiency of the algorithm. Using a standard operator splitting and projection method together with second-order-accurate finite differences  on a uniform three-dimensional mesh, the advancement of the velocity, temperature, and pressure fields at each time step requires the numerical solution of four Helmholtz equations and one Poisson equation. Because these elliptic equations and their boundary conditions are separable, they can be solved efficiently using fast direct methods from the FISHPACK library , with a complexity per problem of $`O(N\mathrm{log}(N))`$, where $`N`$ is the total number of mesh points. Fast direct methods are more efficient than most iterative methods on a single processor , and have the additional advantage that no internal parameters need to be adjusted to obtain convergence. However, fast direct methods are not applicable to complex geometries, to problems with spatially varying parameters, or to complicated boundary conditions that lead to nonseparable equations.

The remainder of this paper is organized as follows. In Section II, we discuss details of our algorithm, namely how the fields and equations are discretized and how the resulting equations are solved. In Section III, we discuss the convergence properties of the algorithm and its efficiency for several representative two- and three-dimensional convection problems. We confirm empirically the second-order accuracy of the solution and examine how the largest time step allowed by stability varies with Prandtl number and with Rayleigh number. Finally, Section IV presents our conclusions and suggests some avenues for further algorithmic improvements. Applications of the algorithm to study quasiperiodic dynamics and spiral defect chaos in three-dimensional boxes can be found in Refs. .

## II Details of the Algorithm

### A Equations and Boundary Conditions

Our goal is to integrate the Boussinesq equations that describe incompressible buoyancy-driven Rayleigh-Bénard convection with an external force $`𝐟`$. These equations can be written in the dimensionless form

$`_tT(t,x,y,z)`$ $`=`$ $`\left[(𝐯\stackrel{}{})T\right]+^2T,`$ (1)
$`_t𝐯(t,x,y,z)`$ $`=`$ $`\left[(𝐯\stackrel{}{})𝐯+\sigma RT\widehat{𝐳}+\sigma 𝐟\right]+\sigma ^2𝐯p,`$ (2)
$`\stackrel{}{}𝐯`$ $`=`$ $`0.`$ (3)
The variables $`x`$ and $`y`$ denote the horizontal coordinates, while the $`z`$ variable denotes the vertical coordinate, with the unit vector $`\widehat{𝐳}`$ pointing in the direction opposite to the gravitational acceleration. The field $`𝐯=(v_x,v_y,v_z)`$ is the velocity field at point $`(x,y,z)`$ at time $`t`$, while $`p`$ and $`T`$ are the pressure and temperature fields respectively. The dimensionless parameters $`\sigma `$ and $`R`$ denote the Prandtl and Rayleigh numbers respectively. The vector field $`𝐟(t,𝐱,T,𝐯)`$ is some external force, e.g., a Coriolis force $`𝐟=2𝐯\times 𝛀`$ arising from a rigid rotation of the convection cell with constant angular velocity $`𝛀=\mathrm{\Omega }\widehat{𝐳}`$. The terms grouped in brackets in Eqs. (1) and (2) are those containing nonlinear terms or linear terms with low-order spatial derivatives, and will be integrated explicitly by the operator splitting method described below.

We would like to integrate Eqs. (1)-(3) in a box geometry defined by the region

$$\frac{\mathrm{\Gamma }_x}{2}x\frac{\mathrm{\Gamma }_x}{2},\frac{\mathrm{\Gamma }_y}{2}y\frac{\mathrm{\Gamma }_y}{2},\frac{1}{2}z\frac{1}{2},$$
(4)
where $`\mathrm{\Gamma }_x`$ and $`\mathrm{\Gamma }_y`$ are the aspect ratios in the $`x`$ and $`y`$ directions respectively (the depth of the fluid has length 1). A no-slip velocity condition on all material walls is assumed

$$𝐯=0,$$
(5)
and the temperature $`T`$ is constant on the bottom and top plates,

$$T=\pm \frac{1}{2}\text{for}z=\frac{1}{2}.$$
(6)
The code allows the temperature on any opposing pair of lateral walls to be periodic, or to satisfy on each lateral wall an arbitrary Dirichlet boundary condition (e.g., a thermally conducting wall corresponding to a linear conducting profile of the form $`T=a+bz`$, where $`a`$ and $`b`$ are constants), or an arbitrary Neumann condition (e.g., a thermally insulating wall with $`_nT=0`$, where $`_n`$ is the normal derivative to the boundary at a given point). To simplify the following discussion, we will consider only the case of insulating sidewalls

$$_nT=0,\text{on lateral walls,}$$
(7)
since the other cases involve just simple modifications. Although the pressure field $`p`$ formally has no associated boundary condition since it does not satisfy a dynamic equation, we will be imposing a Neumann boundary condition on $`p`$ as we explain below (see Eq. (17)).

### B The Time Integration Method

We next discuss the time integration method, since its structure can be explained before having to specify a spatial representation for the fields. In the following subsection, we discuss how the fields and equations are discretized and the latter solved using second-order-accurate finite differences on a uniform spatial mesh.

Our time integration method uses a standard operator splitting and projection method , in which the nonlinear terms containing lower-order or no spatial derivatives are integrated explicitly, then the linear diffusion operators are integrated implicitly, and finally the pressure term $`p`$ is integrated to project the velocity field at the next time step into the space of divergence-free velocity fields. Operator splitting has two benefits. First, the evolution equations for $`T`$ and for the velocity components $`v_i`$ are decoupled from one another, which simplifies the overall algorithm and substantially reduces the total computer memory needed. Second, operator splitting allows larger time steps since the largest time step $`\mathrm{\Delta }t`$ allowed by stability is bounded by a first power of the spatial resolution $`\mathrm{\Delta }x`$, rather than by a second power as would be the case for a fully explicit method.

Let us assume that, at the $`n`$th time step $`t_n=n\mathrm{\Delta }t`$ with $`n0`$, initial fields $`T^n`$ and $`𝐯^n`$ are known that are consistent with the boundary conditions Eqs. (5)–(7). These fields are then advanced to the future values $`T^{n+1}`$ and $`𝐯^{n+1}`$ at time $`t_{n+1}=t_n+\mathrm{\Delta }t`$ as follows:

1. The nonlinear advective term $`N_T[T,𝐯]=(𝐯\stackrel{}{})T`$ of Eq. (1) is integrated explicitly using a second-order-accurate Adams-Bashforth method

$$T^{}=T^n+\frac{\mathrm{\Delta }t}{2}\left(3N_T[T^n,𝐯^n]N_T[T^{n1},𝐯^{n1}]\right),n0,$$
(8)
to produce an intermediate field $`T^{}`$. Here $`T^{n1}`$ and $`𝐯^{n1}`$ denote field values stored from the previous time step $`t_{n1}=t_n\mathrm{\Delta }t`$. For the first time step $`n=0`$ only, a second-order-accurate single-step integrator, Heun’s method , is used in place of Adams-Bashforth to avoid the dependence on the unavailable field values at time $`t=\mathrm{\Delta }t`$.
2. The intermediate field $`T^{}`$ is then advanced to the temperature field $`T^{n+1}`$ at time $`t_{n+1}`$ by using $`T^n`$ as initial data for an implicit Crank-Nicolson step applied to the diffusion term in Eq. (1):

$$\frac{T^{n+1}T^{}}{\mathrm{\Delta }t}=\frac{1}{2}\left(^2T^{n+1}+^2T^n\right).$$
(9)
This can be written as a constant-coefficient Helmholtz equation for the future field $`T^{n+1}`$:

$$\left(1\frac{\mathrm{\Delta }t}{2}^2\right)T^{n+1}=T^{}+\frac{\mathrm{\Delta }t}{2}^2T^n,$$
(10)
and is solved with the boundary conditions Eqs. (6) and (7) applied to $`T^{n+1}`$.
3. Three similar pairs of explicit and implicit steps are then executed successively, first for the velocity component $`v_x`$, then for $`v_y`$, and then for $`v_z`$. The explicit steps advance the velocity field $`𝐯^n`$ to an intermediate field $`𝐯^{}`$, and then the implicit steps advance $`𝐯^{}`$ to a second intermediate field $`𝐯^{}`$. If we denote by $`N_i[T,𝐯]`$ the expressions in brackets of Eq. (2) for $`i=x`$, $`y`$, and $`z`$, then each explicit step has the form

$$v_i^{}=v_i^n+\frac{\mathrm{\Delta }t}{2}\left(3N_i[T,𝐯^n]N_i[T,𝐯^{n1}]\right).$$
(11)
Heun’s single-step method is again used at time $`t_0=0`$ to avoid unavailable field values at time $`t_1=\mathrm{\Delta }t`$. Each field $`v_i^n`$ is next used as initial data for an implicit Crank-Nicolson step that yields a constant-coefficient Helmholtz equation for the field $`v_i^{}`$:

$$\left(1\frac{\sigma \mathrm{\Delta }t}{2}^2\right)v_i^{}=v_i^{}+\frac{\sigma \mathrm{\Delta }t}{2}^2v_i^n.$$
(12)
This is solved with the no-slip boundary condition $`v_i^{}=0`$ on all surfaces, Eq. (5).
4. An incompressible velocity field $`𝐯^{n+1}`$ at time $`t_{n+1}`$ is obtained from the field $`𝐯^{}`$ by integrating the final operator step

$$_t𝐯=p,$$
(13)
with initial data $`𝐯^{}`$, followed by a projection method . We approximate the time derivative in Eq. (13) with a first-order-accurate stencil,

$$\frac{𝐯^{n+1}𝐯^{}}{\mathrm{\Delta }t}=\frac{1}{2}\left(p^{n+1}+p^n\right),$$
(14)
apply the divergence operator to both sides, and then use Eq. (3) in the form $`\stackrel{}{}𝐯^{n+1}=0`$. This yields a Poisson equation for the pressure field $`p`$:

$$^2p^{n+1}=^2p^n+\frac{2}{\mathrm{\Delta }t}\stackrel{}{}𝐯^{}.$$
(15)
Once $`p`$ is known by solving Eq. (15), we obtain $`𝐯^{n+1}`$ from Eq. (14) in the form

$$𝐯^{n+1}=𝐯^{}\frac{\mathrm{\Delta }t}{2}\left(p^{n+1}+p^n\right).$$
(16)
Although mathematically there is no boundary condition for $`p`$—and by discretizing first space and then time, a boundary condition for $`p`$ can be avoided as explained in Refs. —we will solve Eq. (15) with the Neumann condition

$$_np^{n+1}=0\text{on all walls,}$$
(17)
since this allows us to use a fast direct method to solve Eq. (15). There is a substantial literature concerning the appropriateness and accuracy of the boundary condition Eq. (17. Rather than review this literature, we simply point out that a Neumann pressure boundary condition has been shown by previous researchers to produce acceptably accurate results for problems in which the fluid is confined by no-slip surfaces, and we show directly in Section III that our algorithm is second-order accurate in space and second-order accurate in time for several representative problems.

The most time consuming part of this algorithm is, by far, solving the four Helmholtz equations Eq. (10) and Eq. (12) (for $`i=x`$ ,$`y`$, and $`z`$) and solving the Poisson equation Eq. (15).

### C Discretization on a Uniform Mesh

The explicit and implicit steps of the previous section—Eqs. (8) and (10), Eqs. (11) and (12), and Eqs. (15) and (16)—are carried out by discretizing the fields and equations on a single non-staggered mesh of points

$$𝐱_{ijk}=(i\mathrm{\Delta }x,j\mathrm{\Delta }y,k\mathrm{\Delta }z),$$
(18)
that is commensurate with the sides of the box Eq. (4). The mesh indices $`(i,j,k)`$ satisfy

$$\frac{N_x}{2}i\frac{N_x}{2},\frac{N_y}{2}j\frac{N_y}{2},\frac{N_z}{2}k\frac{N_z}{2}.$$
(19)
The aspect ratios $`(\mathrm{\Gamma }_x,\mathrm{\Gamma }_y)`$ and the positive integers $`(N_x,N_y,N_z)`$ are specified as input to the code, and the corresponding spatial resolutions $`(\mathrm{\Delta }x,\mathrm{\Delta }y,\mathrm{\Delta }z)`$ are then determined from the relations $`\mathrm{\Delta }x=\mathrm{\Gamma }_x/N_x`$, $`\mathrm{\Delta }y=\mathrm{\Gamma }_y/N_y`$, and $`\mathrm{\Delta }z=1/(N_z)`$. For large-aspect-ratio convection problems, typically $`\mathrm{\Delta }x=\mathrm{\Delta }y>\mathrm{\Delta }z`$ since the $`x`$ and $`y`$ directions are equivalent and there is a finer structure in the vertical direction caused by the close opposing horizontal plates.

At all mesh points Eq. (18) interior to the box, the first- and second-order spatial derivatives are approximated using centered second-order-accurate 3-point finite-difference stencils. If $`u_{ijk}=u(𝐱_{ijk})`$ denotes the values of a field $`u(𝐱)`$ at the mesh points, then the partial derivative $`_xu`$ at $`𝐱_{ijk}`$ is approximated by

$$\left[_xu\right]_{ijk}\frac{u_{(i+1)jk}u_{(i1)jk}}{2\mathrm{\Delta }x},$$
(20)
with similar expressions for $`_yu`$ and $`_zu`$. The Laplacian of $`u`$ at $`𝐱_{ijk}`$ is approximated by the usual 7-point stencil

$`\left[^2u\right]_{ijk}`$ $``$ $`{\displaystyle \frac{u_{(i+1)jk}2u_{ijk}+u_{(i1)jk}}{\mathrm{\Delta }x^2}}+{\displaystyle \frac{u_{i(j+1)k}2u_{ijk}+u_{i(j1)k}}{\mathrm{\Delta }y^2}}`$ (22)
$`+{\displaystyle \frac{u_{ij(k+1)}2u_{ijk}+u_{ij(k1)}}{\mathrm{\Delta }z^2}}.`$
Nonsymmetric finite-differences are needed to evaluate expressions on those boundaries for which a Neumann condition holds (we will call these “Neumann boundaries”) since field values outside the domain are not available. Thus the right side of the Helmholtz equation, Eq. (10), needs to be evaluated on the Neumann boundaries for which Eq. (7) holds. The value of $`^2T^{}`$ can be approximated there to second-order accuracy by using one-sided 4-point finite-difference approximations for the 2nd-order derivatives, e.g.,

$$[_x^2T^{}]_{0jk}\frac{2T_{0jk}^{}5T_{1jk}^{}+4T_{2jk}^{}T_{3jk}^{}}{\mathrm{\Delta }x^2},$$
(23)
with similar expressions for $`_x^2T^{}`$ at $`x=\mathrm{\Gamma }_x`$, and for $`_y^2T^{}`$ on the $`y=0`$ and $`y=\mathrm{\Gamma }_y`$ boundaries. The divergence $`\stackrel{}{}𝐯^{}`$ on the right side of the pressure equation, Eq. (15), can be approximated to second-order accuracy using the Dirichlet data Eq. (5) and interior field values by replacing Eq. (20) with the following 3-point one-sided finite difference

$$[_xv_x]_{0jk}\frac{3(v_x)_{0jk}+4(v_x)_{1jk}(v_x)_{2jk}}{2\mathrm{\Delta }x}=\frac{4(v_x)_{1jk}(v_x)_{2jk}}{2\mathrm{\Delta }x},$$
(24)
with similar expressions for $`_yv_y`$ and $`_zv_z`$. The advective derivative $`(𝐯\stackrel{}{})T`$ in Eq. (8) vanishes on these Neumann walls since $`𝐯`$ does, and so the explicit steps do not require special treatment.

Given the discretizations Eqs. (20)-(24), the explicit time steps Eqs. (8) and (11) are easily evaluated at all interior points and on the Neumann boundaries. For the implicit steps, the right sides of Eqs. (10) and (12) are also evaluated on the interior mesh points and on the Neumann boundaries. These right sides are then used as input to the FISHPACK  fast direct solvers hw3crt in three dimensions or hwscrt in two dimensions. Also provided as input to the FISHPACK solvers are the corresponding boundary conditions, Eqs. (6) and (7) for $`T`$, Eq. (5) for the velocity components, and Eq. (17) for $`p`$. The FISHPACK solvers return second-order-accurate values (with respect to the spatial resolution) of $`T`$, $`𝐯`$, and $`p`$ on the mesh $`𝐱_{ijk}`$.

We conclude this section with the observation that the discrete velocity field $`𝐯^{n+1}`$ obtained from the concluding step Eq. (16) is only approximately divergence-free even on the mesh points $`𝐱_{ijk}`$, i.e., $`\stackrel{}{}𝐯^{n+1}=O(h^2)`$ where $`h`$ is the larger of the spatial resolutions $`\mathrm{\Delta }x`$, $`\mathrm{\Delta }y`$, and $`\mathrm{\Delta }z`$. This is because the discrete approximation Eq. (20) for the pressure gradient in Eq. (16) is not consistent with the discretization Eq. (22) used to approximate the Laplacian $`\stackrel{}{}`$ on the left side of Eq. (15). The discrete Laplacian can be considered as arising from the evaluation of pressure gradients from pairs of nearest neighbor points as follows

$`[_x^2p]_{ijk}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }x}}\left([_xp]_{(i+1/2)jk}[_xp]_{(i1/2)jk}\right)`$ (25)
$`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }x}}\left({\displaystyle \frac{p_{(i+1)jk}p_{ijk}}{\mathrm{\Delta }x}}{\displaystyle \frac{p_{ijk}p_{(i1)jk}}{\mathrm{\Delta }x}}\right).`$ (26)
In contrast, Eq. (16) evaluates the pressure gradient at a point using a finite difference Eq. (20) that spans three mesh points.

## III Accuracy and Efficiency of the Algorithm

In this section, we discuss several tests that quantify the accuracy of the above algorithm for a convecting fluid in a two-dimensional rectangular domain with periodic sidewalls and in a three-dimensional rectangular domain with perfectly insulating sidewalls, Eq. (7). We first confirm the second-order accuracy of the code with respect to the spatial and time resolutions by studying how the temperature and velocity fields converge with increasing spatial and time resolutions respectively. We next show empirically how the maximum stable time step varies with the Rayleigh number $`R`$ and the Prandtl numbers $`\sigma `$. We then calculate the critical Rayleigh number $`R_c`$ and plot the Nusselt number $`N(R)`$ as a function of the Rayleigh number R, and obtain good agreement with an analytical expression  and with a spectral code . Finally, we show the spatial structure of the fields near onset, to allow comparison with experiment  and with other codes.

We note that, on a workstation with a 667 MHz 21264A 64-bit Alpha processor, a square box with aspect ratio $`\mathrm{\Gamma }=40`$ and spatial resolution $`\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{\Delta }z=1/8`$ takes about 4.8 s per time step of $`\mathrm{\Delta }t=0.001t_v`$. This corresponds to 80 minutes per vertical diffusion time $`t_v`$ and 90 days per horizontal diffusion time $`t_h`$, so this code is too slow to explore $`\mathrm{\Gamma }>20`$ cells over time scales exceeding a horizontal diffusion time. We discuss two ways of improving the efficiency of the code in our concluding comments of Section IV.

### A  Second-Order Convergence With Respect to the Spatial and Time Resolutions

We begin by showing that the order of convergence $`p`$ of the code is asymptotically 2 (second order) in the limits of sufficiently fine spatial and time resolutions. By definition, the convergence with respect to spatial resolution is of order $`p`$ if $`u_hu_{\mathrm{exact}}=O(h^p)`$ in the limit $`h0`$, where $`u=\sqrt{_{ij}u_{ij}^2}`$ denotes the Euclidean norm of a field $`u`$ on the spatial mesh, $`h=\mathrm{\Delta }x=\mathrm{\Delta }z`$ is the uniform spatial resolution in the $`x`$ and $`z`$ directions of a two-dimensional box, $`u_h(x,z)`$ denotes a discrete numerical field on a mesh of resolution $`h`$, and $`u_{\mathrm{exact}}(x,z)`$ is the unknown exact field on the spatial mesh. By writing $`u_h(x,z)=u_{\mathrm{exact}}(x,z)+C(x,z)h^p`$ in the limit $`h0`$, for some function $`C`$ independent of $`h`$, we deduce that the order $`p`$ can be estimated by examining the quantity

$$p_h=\mathrm{log}_2\left(\frac{u_{4h}u_{2h}}{u_{2h}u_h}\right),$$
(27)
in the limit $`h0`$. The estimate Eq. (27) involves field values at the three levels of resolution $`4h`$, $`2h`$, and $`h`$, coarsest to finest. A similar definition for the order of convergence with respect to time resolution can be made if the spatial mesh $`h`$ is replaced by the time step $`g\mathrm{\Delta }t`$,

$$p_g=\mathrm{log}_2\left(\frac{u_{4g}u_{2g}}{u_{2g}u_g}\right),$$
(28)
We first studied the convergence with respect to the spatial mesh $`h`$ for a two-dimensional box with periodic sidewalls, for parameter values $`\mathrm{\Gamma }_x=2\pi /q_c=2.016`$, $`R=1725.01.01R_c`$, and $`\sigma =0.71`$. The initial conditions consisted of a small random perturbation about the linearly conducting state $`T_0=z`$, $`𝐯_0=(u_0,w_0)=\mathrm{𝟎}`$, and these were integrated until a stationary state was attained consisting of two rolls at the critical wave number $`q_c`$. For this small cell, an integration time of $`8t_v`$ was sufficient for the dynamics to become stationary. We then studied the temperature field, $`T_h(x,z)`$, and the $`x`$-component of the velocity field, $`u_h(x,z)`$, for different spatial resolutions $`N=\mathrm{\Gamma }_x/h=16`$, $`32`$, $`64`$, and $`128`$. The time step $`\mathrm{\Delta }t`$ was set respectively to the values $`\mathrm{\Delta }t=0.01`$, $`0.005`$, $`0.0025`$, and $`0.00125`$ since the operating splitting makes the largest stable time step proportional to $`h`$. Table I summarizes the values of the limit Eq. (27) and shows that indeed $`p_h2`$ as $`h0`$, i.e., the code is asymptotically second-order accurate with respect to the spatial resolution $`h`$.

We have also studied the convergence with respect to the time step $`g`$ for a three-dimensional box with perfectly insulating sidewalls and for parameter values $`\mathrm{\Gamma }_x=\mathrm{\Gamma }_y=2`$, $`R=1725.01.01R_c`$, and $`\sigma =0.71`$. The initial condition consisted of small random thermal perturbations, and these were integrated up to 20 diffusion times at which point the state became stationary. For various time resolutions $`g=\mathrm{\Delta }t=0.0001`$, $`0.00005`$, and $`0.000025`$, all with a space resolution of $`N=64`$, the convergence was found (using Eq. (28)) to be $`p=1.68`$. This provides evidence that the code is indeed asymptotically second-order accurate with respect to the time resolution $`g`$.

### B  Dependence of Maximum Stable Time Step on Rayleigh and Prandtl Numbers

Since an important practical feature of any production code is the largest time step that can be taken before numerical instability occurs, we have studied the maximum stable time step as a function of the Rayleigh and Prandtl numbers. A three-dimensional box with periodic sidewalls and aspect ratio $`\mathrm{\Gamma }_x=\mathrm{\Gamma }_y=2`$ was used, with a spatial resolution $`\mathrm{\Delta }x=16`$. The Euclidean norm of the temperature field, $`T`$, was calculated for various values of $`\mathrm{\Delta }t`$ each time after a interval of $`20`$ vertical diffusion times so that transients decayed. The maximum stable time step was then defined as the value of $`\mathrm{\Delta }t`$ such that $`T`$ remains bounded, i.e., $`T<10^5`$.

In Fig. 1(a), we plot the maximum stable time step as a function of Rayleigh number for the two Prandtl number values $`\sigma =1`$ (square symbols) and $`\sigma =10`$ (cross symbols). We see that the maximum stable time step decreases rapidly with increasing Rayleigh number. This is to be expected since the magnitude of the velocity and temperature fields increase with increasing $`R`$. In fact, a best log-log fit to the data yields the relation

$$\mathrm{max}(\mathrm{\Delta }t)R^\alpha $$
(29)
where $`\alpha =1.2`$ when $`\sigma =1`$, and $`\alpha =1.3`$ when $`\sigma =10`$.

In Fig. 1(b), we plot the maximum stable time step as a function of Prandtl number for fixed Rayleigh numbers $`R=2048`$ (square symbols) and $`R=8192`$ (cross symbols). For $`R=2048`$, the maximum stable time step decreases toward both small and large Prandtl numbers. For $`R=8192`$, the maximum stable time step decreases toward small Prandtl numbers but is approximately constant at large Prandtl numbers. The smaller time step needed at small Prandtl numbers can be attributed to the more dynamical nature of the convective flow at small Prandtl numbers, such as the presence of spiral defect chaos .

### C  Estimate of the Critical Rayleigh Number $`R_c`$

A linear stability analysis of the Boussinesq equations about the linearly conducting profile between two infinite horizontal no-slip plates shows that the critical Rayleigh number $`R_c1707.76`$ with critical wave number $`q_c3.117`$, and that the values of $`R_c`$ and $`q_c`$ are independent of the Prandtl number $`\sigma `$ . We tested these predictions and so validated the code by using a two-dimensional box of aspect ratio $`\mathrm{\Gamma }_x=2\pi /q_c=2.016`$, with periodic sidewalls, for Prandtl number $`\sigma =0.71`$. We used a uniform spatial resolution $`h=1/N=\mathrm{\Delta }x=\mathrm{\Delta }z`$ and varied the number $`N`$ of mesh points.

The critical Rayleigh number $`R_c`$ was estimated as the approximate value of R for which the growth rate $`\lambda =\lambda (R)`$ of a small-amplitude (0.01) random perturbation about the linear profile interpolated to zero as a function of R. Thus for a sufficiently tiny initial perturbation of the conducting profile, there is a time interval over which the $`z`$-velocity component $`w`$ grows approximately exponentially

$$w(t,x,z)c(R)e^{\lambda t},$$
(30)
where $`\lambda `$ is the growth rate, and $`c`$ is independent of $`t`$ but can vary with R. For $`R>R_c`$, the growth rate is positive, for $`R<R_c`$, the growth rate is negative and interpolating between known positive and negative values provides an estimate of $`R_c`$, for which $`\lambda =0`$.

Our protocol was to set $`R=R_+=1725=1.01R_c`$ just above onset, set the initial velocity field to zero, $`𝐯_0=(u_0,w_0)=\mathrm{𝟎}`$, and set the initial temperature field $`T_0=z+\delta T`$ to a tiny random perturbation $`\delta T(x,z)`$ of the linear profile $`T=z`$, with $`|\delta T|0.01`$. The initial conditions were then integrated for a short time and the growth rate estimated from the formula

$$\lambda _+\frac{\mathrm{ln}\left(w(t_2,x,z)/w(t_1,x,z)\right)}{t_2t_1},$$
(31)
where $`t_2`$ and $`t_1<t_2`$ are two times during the exponential growth of the magnitude of the $`z`$-component of the velocity field $`w`$. The calculation was then repeated with the same initial condition but for $`R=R_{}=1691=0.99R_c`$ to estimate a decay rate $`\lambda _{}`$. The critical Rayleigh number was then estimated as the zero of the line joining the points $`(R_+,\lambda _+)`$ and $`(R_{},\lambda _{})`$. The estimated critical Rayleigh numbers $`R_c`$ as a function of the number of mesh points $`N`$ are summarized in Table II. The values are correct to a relative error of better than one percent for the finest spatial resolution, confirming the correctness and convergence of the discretization and of the solution technique.

### D The Nusselt number versus Rayleigh number Curve $`N(R)`$

Another way to characterize the accuracy of a convection code is by the dimensionless Nusselt number $`N(t,R,\sigma )`$, which is the instantaneous global vertical heat transport through the fluid layer, normalized to the heat transport arising from thermal conduction alone. For the dimensionless variables used in Eqs. (1)-(3) above, the Nusselt number can be expressed in the form

$$N=1+w(TT_{\mathrm{cond}}),$$
(32)
where $`w`$ is the $`z`$-component of the velocity field and $`T_{\mathrm{cond}}=z`$ is the temperature profile of the linear conducting state with $`𝐯=\mathrm{𝟎}`$. The brackets $`\mathrm{}`$ denote an average of a quantity over the horizontal coordinates. Sufficiently close to onset, numerical values of $`N`$ can be compared with an analytical expression  that is valid asymptotically in the limit $`RR_c0^+`$.

We have evaluated Eq. (32) for the two-dimensional domain with periodic sidewalls of the previous section, with parameters $`\mathrm{\Gamma }_x=2.016`$, $`\sigma =0.71`$, and $`N=16`$. Starting from a small perturbation of the linear profile, we integrated until a stationary state was attained, and then evaluated Eq. (32) for the stationary state. Fig. 2 shows how $`N`$ empirically varies with R, and we compare this curve with the analytical result of Schlüter et al , and with numerical values obtained by Clever and Busse, who used a two-dimensional spectral code with periodic side walls. The agreement is good in both cases and confirms the correctness and accuracy of the code.

### E Spatial Structure of the Numerical Solutions

We conclude this section with a few examples of the spatial structure of the fields obtained from the code, to show that the fields are physically reasonable when adequately resolved, and are qualitatively in agreement with other codes and with experiment .

for $`\sigma =0.71`$, Fig. 3 shows contours of constant temperature $`T`$ and for two values of the Rayleigh number, $`R=2500`$ in (a) and $`R=10^4`$ in (b). Warm fluid ascends in the middle of the cell and descends as cooler fluid on both sides. As $`R`$ increases, a thermal boundary layer forms at the top and bottom plates, creating a finer spatial structure that will require eventually a decrease in the vertical spatial mesh size $`\mathrm{\Delta }z`$. Fig. 4(a) shows the corresponding velocity field $`𝐯=(u,w)`$, while Fig. 4(b) shows the vertical component $`w`$ through the midline of the cell. The occurrence of two square-shaped convection cells of opposite vorticity is in good agreement with experiment .

Fig. 5 shows constant temperature contours in a three-dimensional box with insulating sidewalls at time $`t=200t_v`$, for parameters $`\mathrm{\Gamma }=16`$, $`R=2500`$, $`\sigma =0.71`$, and $`h=\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{\Delta }z=1/8`$, $`\mathrm{\Delta }t=0.01`$. In agreement with experiment  and with calculations on the Swift-Hohenberg model of convection , the rolls are approximately normal to the lateral walls and the pattern consists of two diagonally opposite foci. For slightly higher $`R=8500`$, Fig. 5(b) shows that the oscillatory instability commenced in the form of ripples that propagate along the length of the rolls. The occurrence of the oscillatory instability and its spatial form are in good agreement with the linear stability analysis of Busse and collaborators  and with experiment.

## IV Conclusions

We have described and characterized a semi-implicit finite-difference algorithm for integrating the Boussinesq equations in two- and three-dimensional boxes, with sidewalls that are periodic, thermally insulated, or thermally conducting. Our approach is useful for simple geometries like a box, cylinder, torus, and annulus, with boundary conditions such that various linear operators are separable so that fast direct methods can be applied. The resulting algorithm is sufficiently efficient that aspect ratios up to $`\mathrm{\Gamma }20`$ can be studied on a single-processor workstation over several days. We verified that the code was second-order-accurate with respect to the spatial and time resolutions, and that it gave good agreement for the critical Rayleigh number and for the Nusselt number versus Rayleigh number curve near onset.

The most significant feature of our algorithm is the use of a single non-staggered mesh for discretizing the equations and fields (velocity, temperature, and pressure). The use of a single mesh simplifies the writing and validation of the code, and facilitates adding new physical terms like a Coriolis force. The single mesh also allowed the use of fast direct methods from the FISHPACK library  to solve the Helmholtz and Poisson equations associated with the implicit part of each time step. We found that numerical integrations of the Boussinesq equations were stable on a single mesh despite results of some previous papers that suggested that a non-staggered Navier-Stokes code could be unstable because of pressure oscillations,

Although the algorithm is useful and has been successfully applied to several problems , there are two ways that the algorithm could be improved for the future study of large-aspect-ratio Rayleigh-Bénard convection. First is to parallelize the code for a distributed memory parallel computer so that aspect ratios comparable to the largest experiments ($`50<\mathrm{\Gamma }<100`$) could be studied. This is technically straightforward and would involve, first, distributing the arrays that represent the fields over the various processors, and, second, replacing the fast direct solvers with iterative methods for sparse matrices. Because of the simpler data structures and reduced communication overhead associated with a finite-difference discretization, the parallelized algorithm will likely be more efficient for simple geometries and for simple boundary conditions than the parallel spectral element method of Refs. .

Because time integration algorithms involve sequential steps, parallelizing a code allows a larger spatial domain, but not a longer observation time, to be studied for a fixed amount of wall-clock time. A second helpful improvement would be to increase the efficiency of the time integration method close to the onset of convection so that larger time steps can be taken for a given computational effort. A weakness of the operating splitting used in most convection codes—finite difference, spectral, and spectral-element—is that the explicit integration of the advection terms imposes a bound on the time step of the form $`Cϵ^{1/2}\mathrm{\Delta }x`$ where $`C`$ is a constant and $`ϵ=(RR_c)/R_c`$ is the reduced Rayleigh number. This bound is independent of the spatial resolution and diverges less rapidly in the limit $`ϵ0^+`$ than the physical time scale, which is proportional to $`ϵ^1`$ . It would be interesting to explore whether a more sophisticated explicit time-stepping technique such as a matrix exponential method  or a fully implicit method  may succeed in allowing larger time steps that are commensurate with the physical time scale while retaining the efficiency of the present code.

## Acknowledgments

We would like to thank Michael Cross, Paul Fischer, and Mark Paul for helpful discussions, and the Department of Energy for supporting this research under grant DE-FT02-98ER14892.

Figure 1.

Figure 2.

Figure 3.

Figures 4.

Figure 5.