Low-energy local density of states of the 1D Hubbard model

Interacting one-dimensional quantum wires are well-studied examples of systems in which the Fermi liquid paradigm of electron-like quasi-particles is known to break down. Luttinger liquid theory predicts that strong correlations lead to the remarkable phenomenon of separate spin- and charge-density waves as the fundamental collective excitations in one dimension at low energies [1]. Experimental confirmation for this picture has long been controversial, but now there is some evidence for separately dispersing spin and charge resonances for quantum wires on semiconductor hetero-structures [2, 3], quasi–one-dimensional crystals [4] and self-organized atomic chains [5] as a function of momentum. The density of states (DOS) has also been analyzed by scanning tunneling spectroscopy (STS) in carbon nanotubes [6–8] and self-organized atomic chains [9]. This immediately invites the question as to whether there are characteristic signatures from standing waves of separate spin and charge densities, that can potentially be detected in locally resolved STS experiments in finite wires. Detailed Luttinger liquid calculations near boundaries exist, which predict corresponding wave-like modulations in the local DOS and a sharp reduction of the DOS near boundaries [10–16]. However, it is far from clear whether these signatures are robust in realistic lattice systems, since even minimalistic models like the Hubbard chain are not perfect Luttinger liquids. There are two main reasons for possible discrepancies: First, the assumed degeneracy of spin and charge modes in the low-energy spectrum can never be exact and will be lifted by band curvature and other effects. Second, it is known that strong logarithmic corrections will arise from a spin umklapp operator which is generically present in 1D electron systems with SU(2) invariant interactions [17–19].

The most relevant minimalistic lattice model for quantum wires is the Hubbard chain,

$$\begin{equation} H = -t \sum_{\sigma,\, x=1}^{L-1} \left( \psi_{\sigma, x}^\dagger \psi_{\sigma, x+1}^{\phantom\dagger} + {\rm h.c.} \right) + U \sum_{x=1}^L n_{\uparrow, x} n_{\downarrow, x}, \end{equation} \tag{ 1 }$$

which captures the main aspects of interacting one-dimensional electron systems. This model shows signatures of spin-charge separation in numerical simulations of the momentum resolved DOS [20], which is the central quantity relevant for photoemission experiments. In this paper we will now analyze the local DOS as a function of position and energy which in turn is relevant for STS experiments.

Recently, several numerical density matrix renormalization group (DMRG) studies considered the local DOS for spinless lattice fermion models in one dimension [21–23], where there is no spin and charge separation. Therefore, the complications mentioned above do not arise [24, 25] and the agreement with theory is close to perfect in that case. As we will show here, the local DOS for a spinful electron system away from half-filling is much more complex and some key features of the Luttinger liquid theory are strongly renormalized.

The observable of interest is the local DOS to tunnel an electron with spin up at a certain energy ω into the wire,

$$\begin{equation} \rho(x, \omega) = \sum_\alpha \, \left| \langle \left. \omega_\alpha \right| \psi_{\uparrow,x}^\dagger \left| 0 \right. \rangle \right|^2 \,\delta(\omega - \omega_\alpha) \end{equation} \tag{ 2 }$$

$$\begin{equation} = \frac{1}{\pi}{\rm Im} \int_0^\infty e^{\, \mathrm i \omega t} \, \mathrm i \langle \psi_{\uparrow,x}^{\phantom\dagger}(t) \psi_{\uparrow,x}^\dagger(0) \rangle \, \mathrm d t . \end{equation} \tag{ 3 }$$

Here, the fermion operator ψ^†_↑,x creates a particle at position x on top of the ground state $\left | 0 \right . \rangle$ and the sum runs over all states $| \omega _\alpha \rangle$ with one additional electron.

In the following we will compare analytic calculations for the local DOS from Luttinger liquid theory with simulations using the numerical DMRG algorithm [26, 27]. For the matrix elements in eq. (2) we use a multi-target DMRG with a large number of target states in two different particle number sectors [21] and keep track of all local fermion creation operators.

In order to identify the separate spin and charge excitations, we first focus on the energy spectrum Δω = ω_α − ω₀ for finite wires of length L as a function of interaction U as shown in fig. 1 in units of t = 1. The ground state |0〉 is assumed to be a filled Fermi sea with total magnetization S^z = 0, so that all particle excited states |ω_α〉 have N_↑ = N_↓ + 1 particles and total spin z-component S^z = 1/2.

Zoom In Zoom Out Reset image size

For U = 0 all excitations are described by simple products of fermion operators $c_{\sigma ,n}^\dagger\kern-1pt =\kern-1pt \sqrt {\kern-1pt\frac {2}{L+1}}\kern-1pt\sum _x\kern-1pt \psi _{\sigma ,x}^\dagger \sin (k_F\kern-1.5pt+\kern-1.5pt k_{n+1})x$, where $k_n= n\frac {\pi }{L+1}$ and |ω₀〉 = c^†_↑,0|0〉 is the lowest-energy particle excitation. Since the spectrum is approximately linear Δω ∼ v_F (k − k_F ), a state involving several fermion operators (i.e., a multi-particle excitation) is nearly degenerate with a single-particle excitation at U = 0 if the sum of excited wave numbers is the same, which results in the quantized fermion levels shown in fig. 1. For example, in level 1, the two states c^†_↑,1|0〉 and c^†_↑,0c^†_↓,0c _↓,−1|0〉 both have approximately the same energy ω₁ − ω₀ ≈ v_F k₁. In general there are many multi-particle states in each level n, but only one single-particle state c^†_↑,n|0〉 carries all the DOS if U = 0.

The situation changes for finite U. Now all states may potentially carry spectral weight in the DOS and the near-degeneracy of the levels is lifted. According to the Luttinger liquid picture the states are now described by integer spin and charge quantum numbers {m_s,m_c} with energies $\omega _{m_s,m_c}=(m_{s}v_s+m_c v_c) \frac {\pi }{L+1}$ in terms of the spin and charge velocities v_s ⩽ v_c [1, 14]. The spectrum in fig. 1 shows a regular spin and charge spacing with quantum numbers shown on the right. For each spin/charge mode {m_s,m_c} there can be several states with approximately the same energy. The number of states in each mode is given by the number of ways it can be created by boson creation operators b^†_m,ν, i.e., the product of the integer partitions of m_s and m_c. For example, for the mode 0s 2c in fig. 1, we have m_c = 2 = 1 + 1 corresponding to the two states created by b^†_2,c and (b^†_1,c)², respectively, which indeed have almost the same energy. On the other hand, the near-degeneracies of spin modes (e.g., 2s 0c) are significantly split in fig. 1 by well-known logarithmic corrections as will be discussed below [17–19].

Regarding the DOS, it is now interesting to explore how the total spectral weight is distributed among the excited states. The compact answer from Luttinger liquid theory is that the summed up DOS in each mode {m_s,m_c} should be proportional to a power law of the corresponding energy, but there are no predictions how the DOS is distributed among the states within each mode. For the lowest-energy modes the situation is still simple, since the first level corresponds to one spin and one charge state, which should have roughly the same DOS of 1/2 each for small U according to theory. However, the numerical results in fig. 2 already demonstrate obvious deviations from this prediction, since the charge mode has a much larger DOS for small U. The reason for this discrepancy comes from the non-linear band curvature, which lifts the degeneracy for finite L even at U = 0, so that the interaction has to overcome this energy splitting. Indeed for cases in which the two states are exactly degenerate at U = 0 (e.g., in the thermodynamic limit) the spin and charge states have comparable DOS even for infinitesimally small U. As can be seen in fig. 2 for L = 90 the charge state dominates for U ≲ 0.5, which would imply that for $U\lesssim 10 \frac {v_s \pi }{L}$ the band curvature is dominant over the interaction effects. The next five states from the second fermion level in fig. 2 show a similar crossover behavior with U. Nonetheless, features in the local DOS will clearly show the interaction effects even in the crossover region as we will see below.

Zoom In Zoom Out Reset image size

In the 0s 2c and the 2s 0c modes there are two states each as expected. However, only one of the states in each mode contributes most of the corresponding DOS, while the other state would be practically invisible in an STS experiment. The reason why some states in a given mode have zero DOS can sometimes be linked to exact symmetries, such as the SU(2) symmetry of generic Coulomb interactions. In particular, all particle excitations |ω_α〉 in the S^z = 1/2 sector are representatives of SU(2) multiplets. The lowest-energy state |ω₀〉 always belongs to a doublet with S = 1/2. Excitations with charge bosons b^†_ℓ,c never change the total spin, but excitations with spin bosons b^†_ℓ,s may generate higher spin values, which can be calculated using the commutation rules of the non-Abelian SU(2)-Kac-Moody algebra [17], since the spin bosons correspond to the modes of the SU(2) current along the z-direction as discussed in the appendix. For example the lowest-energy S = 3/2 state is given by the spin boson excitation $\frac {1}{\sqrt {3}}[(b^\dagger _{1,s})^2 - b^\dagger _{2,s}] |\omega _0\rangle$ with m_s = 2 which is plotted as a dashed line in fig. 1. The total spin analysis is discussed in more detail in the appendix. It is useful since states with S > 1/2 must have zero DOS due to angular momentum addition rules. For example in the fifth spin mode 5s 0c, there are 7 states, three of which have total spin S = 3/2 and exactly zero DOS. Interestingly, three more states carry only very small spectral weight, so that only one state dominates for this mode. This indicates that the eigenstates remain the same in terms of their bosonic expressions even for U ≠ 0. For charge modes there is also exactly one state which carries the overwhelming weight in each mode, but the DOS of the other states is finite and generally increases with U.

For a complete analysis it is now useful to turn to the local DOS ρ_ms,mc(x) for each mode {m_s,m_c}. For a non-interacting system, the local DOS is given by the square of a quickly oscillating standing wave, corresponding to the eigenstate at the corresponding wave vector near k_F . However, with interactions we also expect long-wavelength modulations and a characteristic behavior near the boundary. It is therefore useful to decompose the local DOS into a quickly oscillating part (o) and a uniform part (u). In order to analyze the long-wavelength modulations we will concentrate on the uniform part of the local DOS in the following, which according to bosonization is given by a product of spin and charge contributions [14]

$$\begin{equation} \rho_{m_s,m_c}(x) = |c_x|^2 \rho^{\rm u}_{s,m_s}(x)\rho^{\rm u}_{c,m_c}(x) \end{equation} \tag{ 4 }$$

for a given mode {m_s,m_c}. The slowly varying amplitudes ρ^u_ν,m(x) can be determined by a simple recursive formula for spin and charge (ν = c,s) separately [14],

$$\begin{equation} \rho^{\rm u}_{\nu, m}(x) = \frac{1}{m} \sum_{\ell=1}^m \rho^{\rm u}_{\nu, m-\ell}(x) \gamma^{\rm u}_{\nu,\ell}(x), \end{equation}\noindent \tag{ 5 }$$

where

$$\begin{equation} \gamma^{\rm u}_{\nu,\ell}(x)=a_\nu + b_\nu \, \cos(2k_\ell x). \end{equation} \tag{ 6 }$$

Here we have defined spin and charge exponents a_ν = (1/K_ν + K_ν)/4 and b_ν = (1/K_ν − K_ν)/4 in terms of the respective Luttinger parameters K_ν for ν = c,s. The overall prefactor $|c_x|^2\propto (\sin \frac {\pi x}{L+1})^{b_c+b_s}$ in eq. (4) does not depend on energy and serves as normalization so that ρ^u_ν,m=0 = 1. It is straightforward to see that the recursive formula results in power laws for the DOS in the bulk ρ∝ω^a_c+a_s−1 for L → ∞ [14, 28]. In addition, the formula predicts slow wave-like modulations in the local DOS due to the second term in eq. (6), which also survive in the thermodynamic limit near the edge [10].

The local DOS from the DMRG data is shown in fig. 3 for the first few modes at 1/3-filling. For the lowest excitation |ω₀〉 the oscillating and uniform parts are the same and given by the prefactor |c(x)|². Already at first sight it is surprising to see that the local DOS for all modes increases slightly near the boundary, while all previous analytic calculations have predicted it to decrease according to the boundary exponent [10, 28, 29]. Indeed it must be emphasized that the local DOS does not fit the theoretical prediction. All curves should in principle be fit free, up to one overall normalization, since the local DOS of all levels follows from eqs. (4)–(6), where the Luttinger parameters K_c(U) and K_s = 1 are known from the thermodynamic Bethe ansatz. However, in fig. 3 two important adjustments have been made: First the theoretical curves were shifted down for the charge modes and up for the spin modes in order to fit the numerical data (indicated by arrows). This adjustment was already observed in the crossover of the total DOS in fig. 2 due to the competition of energy scales (band curvature vs. interaction) as argued above. We could not find any explanation of the shifts in the framework of a renormalized Luttinger liquid. Secondly, we find that the spin Luttinger liquid parameter must be chosen considerably larger than unity K_s ≈ 1.16 for all spin modes in order to fit the numerical data corresponding to attractive behavior in the spin modes. This is especially surprising since the charge Luttinger parameter from the Bethe ansatz K_c = 0.9081 agrees perfectly with the data without any finite-size adjustments. There are no other adjustable parameters in the fits of fig. 3, except for one overall normalization constant. The oscillating parts ρ^o_ν can be analyzed analogously and give the same results (not shown). The local DOS can also be calculated by Fourier transforming the real-space Green's function [10–13, 28], but for an analytic analysis of individual states the recursive formulas in eqs. (4)–(6) are a great advantage, since they give a closed form without divergences and in return allow the quick and accurate determination of the effective parameters K_s and K_c by fitting the data from (numerical) experiments for the first spin and charge excited states.

Zoom In Zoom Out Reset image size

In fig. 4 the behavior of the Luttinger parameters from the corresponding fits to the uniform local DOS is shown as a function of interaction and length at 1/3-filling. The charge parameter from the Bethe ansatz K_c always agrees very well with the data without any additional adjustments. However, the observed spin parameter K_s is considerably larger than K_s = 1. Non-Abelian bosonization predicts K_s = 1 for any SU(2) invariant model, but at the same time it is known that a marginal irrelevant operator causes corrections to the anomalous dimension which only vanish logarithmically slowly with 1/ln L in the thermodynamic limit [17, 18]. In Abelian bosonization such a correction can indeed effectively be modeled by a renormalizing spin Luttinger parameter [17–19, 30, 31]

$$\begin{equation} K_s-1\propto g, \quad {\rm with }\quad g^{-1} + \frac{1}{2} \ln(g) = \ln(L/L_0), \end{equation} \tag{ 7 }$$

where L₀ is non-universal and depends on the model and the quantity of interest. As shown in fig. 4 such a renormalization description is indeed consistent with our data for K_s. The parameter K_s increases with U at a given length L, but decreases slowly as the length is increased. The parameter K_s appears to be the same for all spin modes at a given U and L, i.e., independent of energy ω. The renormalization of K_s is very slow, so that exponentially large systems are required to observe the thermodynamic limit K_s → 1. The fit parameters for K_s ≈ 1 + 1.98g with ln L₀ ≈ −6 are outside the range, which would normally be expected for a spin chain model [19]. Therefore, the particular form of the observed corrections remains a puzzle.

Zoom In Zoom Out Reset image size

Nonetheless, the results of the logarithmic corrections have interesting consequences. In particular, the corrections are so large, that the boundary exponent α_B = (1/K_s + 1/K_c)/2 − 1 may become negative if K_s + K_c < 2K_sK_c, i.e., K_s − 1 ≳ 1 − K_c to lowest order in the correction, which is indeed the case for small U ≲ 2t at L = 92. This results in an increase of the local DOS for small energies ω and small distances x near the boundary, which is described by a weak power law divergence ρ∝x^b_s+b_cω^α_B. Such a negative boundary exponent would also explain the recently reported boundary anomalies for small interactions in functional renormalization group studies [32, 33]. It is known that logarithmic corrections also exist as a function of energy near the boundary, which lead to a very small energy region where an uncorrected boundary exponent can be observed [24, 25].

In conclusion, we have analyzed the local DOS of the Hubbard model in the low-energy regime. Individual states can be classified by separate spin and charge quantum numbers. We observe that typically only one eigenstate has a relatively large local DOS in each spin/charge mode, while all other states in that mode are negligible. The spin and charge Luttinger parameters K_s and K_c can be extracted from the modulations in the local DOS of individual excited states. While the charge parameter K_c agrees well with the Bethe ansatz, the spin Luttinger liquid parameter is attractive K_s > 1 due to large finite-size corrections, which can only be neglected for exponentially large chain lengths. In fact, the corrections to K_s are unexpectedly strong and may even lead to negative boundary exponents for moderate interactions U ≲ 2t. The common assumption that it is possible to generically use K_s = 1 due to SU(2) invariance is certainly not justified for the local DOS. In particular, for finite wires on conducting substrates the interactions may be reduced by screening, so that the charge Luttinger liquid parameter K_c may be close to unity, while the spin Luttinger liquid parameters K_s can already be significantly increased, which leads to a negative boundary exponent. This would have quite dramatic consequences, since the DOS near the boundary determines the tunneling between connected wires and the renormalization of the conductivity through impurities [29], which will show an increase at low temperatures in this scenario.

We are thankful for useful discussions with A. Struck and M. Bortz. This work was supported by the DFG and the State of Rheinland-Pfalz via the SFB/Transregio 49 and the MAINZ graduate school of excellence.

In order to determine the total spin eigenstates as shown in tables 1 and 2, it is useful to use non-Abelian bosonization in the spin channel [17]. In this case the excitations are created by the modes of SU(2) currents J^a_m with a = x,y,z obeying the Kac-Moody algebra,

$$\begin{equation} [{J^a_m},{J^b_n}] = \, \mathrm i \varepsilon^{abc} J^c_{m+n} + \textstyle\frac{1}{2}m \delta^{a,b} \delta_{m,-n}. \end{equation} \tag{ A.1 }$$

The ground state is characterized by J^a_m|ω₀〉 = 0, ∀ m < 0. The total spin operator is given in terms of the m = 0 currents,

$$\begin{equation} S^2 = \vec{J}_0 \cdot \vec{J}_0 = 2J_0^+J_0^- + J_0^z + (J_0^z)^2, \end{equation} \tag{ A.2 }$$

where J^± = J^x ± iJ^y. The current modes in the z-direction are related to the Abelian spin bosons above by $J^z_m = \sqrt {\frac {m}{2}} b^\dagger _{m,s}$ and $J^z_{-m} = \sqrt {\frac {m}{2}} b_{m,s}$ for m > 0. It is therefore straightforward to consider the total spin of any bosonic spin and charge excitation by using the Kac-Moody commutation relations. The charge bosons commute with the total spin operator S². Spin excitations are created by products of spin creation operators b^†_m,s acting on |ω₀〉. The corresponding normalized states can be labelled by the set of which bosons were created |{m₁,m₂,m₃,...}〉, e.g., $|\{3,1,1\}\rangle = \frac {1}{\sqrt {2}} b^\dagger _{3,s} (b^\dagger _{1,s})^2 |\omega _0\rangle$. Therefore, the matrix elements of 〈{m₁,m₂,m₃,...}|S²|{m₁',m₂',m₃',...}〉 between any two such excitations can be evaluated uniquely by the Kac Moody algebra (A.1). The J^z₀ operators commute with all excitations and the ground state is characterized by $J^z_0 |\omega _0\rangle = S^z|\omega _0\rangle = \frac {1}{2}|\omega _0\rangle$ in our case. For the J^± operators we use the Kac-Moody relation in eq. (A.1) with the help of computer algebra in order to successively commute them to the right until the action on the ground state is known. This results in a non-diagonal matrix for S² for each spin mode separately, which can be brought into diagonal form. The resulting eigenstates and eigenvalues are given in table 1 for the case in which the state |ω₀〉 has total spin of S = 0. In this paper we considered an excitation on the filled Fermi sea, which will correspond to a state |ω₀〉 with S = 1/2. In this case the eigenstates are given in table 2.

Table 1:. The SU(2) eigenstates in terms of bosonic quantum numbers with integer spin.

S=0	${\left | \omega _0 \right\rangle }$
S=1	$\left | \{1\} \right\rangle$
S=0	$\left | \{1,1\} \right\rangle$
S=1	$\left | \{2\} \right\rangle$
S=0	$\left | \{2,1\} \right\rangle$
S=1	$\left | \{3\} \right\rangle$
S=1	$\left | \{1,1,1\} \right\rangle$
S=0	$\sqrt {\frac {2}{3}} \left | \{3,1\} \right\rangle + \sqrt {\frac {1}{3}} \left | \{1,1,1,1\} \right\rangle$
S=0	$\sqrt {\frac {1}{9}} \left | \{3,1\} \right\rangle + \sqrt {\frac {2}{3}} \left | \{2,2\} \right\rangle - \sqrt {\frac {2}{9}} \left | \{1,1,1,1\} \right\rangle$
S=1	$\left | \{4\} \right\rangle$
S=1	$\left | \{2,1,1\} \right\rangle$
S=2	$-\sqrt {\frac {2}{9}} \left | \{3,1\} \right\rangle + \sqrt {\frac {1}{3}} \left | \{2,2\} \right\rangle + \sqrt {\frac {4}{9}} \left | \{1,1,1,1\} \right\rangle$
S=0	$\sqrt {\frac {4}{7}} \left | \{4,1\} \right\rangle + \sqrt {\frac {3}{7}} \left | \{2,1,1,1\} \right\rangle$
S=0	$\sqrt {\frac {6}{63}} \left | \{4,1\} \right\rangle + \sqrt {\frac {49}{63}} \left | \{3,2\} \right\rangle - \sqrt {\frac {8}{63}} \left | \{2,1,1,1\} \right\rangle$
S=1	$\left | \{5\} \right\rangle$
S=1	$\left | \{3,1,1\} \right\rangle$
S=1	$\left | \{2,2,1\} \right\rangle$
S=1	$\left | \{1,1,1,1,1\} \right\rangle$
S=2	$-\sqrt {\frac {1}{3}} \left | \{4,1\} \right\rangle + \sqrt {\frac {2}{9}} \left | \{3,2\} \right\rangle + \sqrt {\frac {4}{9}} \left | \{2,1,1,1\} \right\rangle$

Table 2:. The SU(2) eigenstates in terms of bosonic quantum numbers with half-integer spin.

S=1/2	${\left | \omega _0 \right\rangle }$
S=1/2	$\left | \{1\} \right\rangle$
S=1/2	$\sqrt {\frac {2}{3}} \left | \{2\} \right\rangle + \sqrt {\frac {1}{3}} \left | \{1,1\} \right\rangle$
S=3/2	$-\sqrt {\frac {1}{3}} \left | \{2\} \right\rangle + \sqrt {\frac {2}{3}} \left | \{1,1\} \right\rangle$
S=1/2	$\sqrt {\frac {1}{3}} \left | \{3\} \right\rangle + \sqrt {\frac {2}{3}} \left | \{1,1,1\} \right\rangle$
S=1/2	$\sqrt {\frac {2}{9}} \left | \{3\} \right\rangle + \sqrt {\frac {6}{9}} \left | \{2,1\} \right\rangle - \sqrt {\frac {1}{9}} \left | \{1,1,1\} \right\rangle$
S=3/2	$-\sqrt {\frac {4}{9}} \left | \{3\} \right\rangle + \sqrt {\frac {3}{9}} \left | \{2,1\} \right\rangle + \sqrt {\frac {2}{9}} \left | \{1,1,1\} \right\rangle$
S=1/2	$-\sqrt {\frac {1}{3}} \left | \{3,1\} \right\rangle + \sqrt {\frac {2}{3}} \left | \{1,1,1,1\} \right\rangle$
S=1/2	$\begin{array}{@{}l@{}}\sqrt{\frac {12}{27}} \left | \{4\} \right\rangle + \sqrt {\frac {2}{27}} \left | \{3,1\} \right\rangle + \sqrt {\frac {12}{27}} \left | \{2,1,1\} \right\rangle\\ \quad+\sqrt{\frac{1}{27}} \left | \{1,1,1,1\} \right\rangle\end{array}$
S=1/2	$\begin{array}{@{}l@{}}\sqrt{\frac {3}{54}} \left | \{4\} \right\rangle + \sqrt {\frac {8}{54}} \left | \{3,1\} \right\rangle + \sqrt {\frac {27}{54}} \left | \{2,2\} \right\rangle\\ \quad- \sqrt {\frac {12}{54}} \left | \{2,1,1\} \right\rangle + \sqrt {\frac {4}{54}} \left | \{1,1,1,1\} \right\rangle\end{array}$
S=3/2	$-\sqrt {\frac {1}{3}} \left | \{4\} \right\rangle + \sqrt {\frac {1}{3}} \left | \{2,2\} \right\rangle + \sqrt {\frac {1}{3}} \left | \{2,1,1\} \right\rangle$
S=3/2	$\begin{array}{@{}l@{}}-\sqrt {\frac {1}{6}} \left | \{4\} \right\rangle + \sqrt {\frac {4}{9}} \left | \{3,1\} \right\rangle - \sqrt {\frac {1}{6}} \left | \{2,2\} \right\rangle\\ \quad + \sqrt {\frac {2}{9}} \left | \{1,1,1,1\} \right\rangle\end{array}$