\section{Introduction}
Loop quantum gravity is an background independent way to quantize gravity, where the Hamiltonian constraint  is written in terms of
Ashtekar-Barbero connection and densitized triad \cite{1}-\cite{1a}.
In loop quantum gravity the curvature of the connection is
expressed through
the holonomy around a loop. The area of such a loop cannot take
arbitrary small nonzero values because
 the area operator in loop quantum gravity has a discrete spectrum.
 At a kinematics level
the Hilbert space of loop quantum gravity is a space of spin networks
\cite{s}-\cite{sa}, and the time evolution of these spin
networks naturally
leads to the spin foam model \cite{f}-\cite{fa}.
 Even though in loop quantum gravity the geometry is dynamical, the topology is fixed.
This is because, like the canonical formalism, loop quantum gravity is a second quantized formalism.
Just like  the processes where the particle number is not conserved
cannot be fully analysed in the first quantization formalism, and we have to resort to second quantization to analyse such processes;
 the topology changing processes  cannot be fully analysed in second quantized formalism,
 and we have to resort to third quantization  to analyse such processes
 \cite{t}-\cite{ta}.
A third quantization of loop quantum gravity naturally
leads to group fields theory,
 which is a field theory defined on a group manifold \cite{g}-\cite{ga}.
 Topology changing process can be analysed in the framework of group field theory. These processes occur due to the creation and annihilation
of spin states. So, the Feynman amplitudes of group field theory can be
equivalently expressed in terms of spin-foam models.
Many models for  loop quantum cosmology  have
 been recently analysed and are
 better understood than full loop quantum gravity \cite{c}-\cite{ca}.
 Similarly, group field cosmology has been studied
 as a mini-superspace approximation to the
full group field theory \cite{m}-\cite{ma}. In group field cosmology first
the Hamiltonian constraint of loop quantum cosmology
is viewed as the kinetic part of the group field theory and then  interactions are added.
The bosonic group field cosmology has already been studied and
 there is no reason to restrict the statistics of the group field cosmology to bosonic case. Hence, in this paper we will also analyse a fermionic
group field cosmology. We will in fact construct a supersymmetric group field cosmology in ${\mathcal{N}}=1$ superspace formalism with gauge symmetry.
Even though  color has already been
incorporated in group field theory \cite{color}-\cite{colora}, it has not been so far been incorporated in any model of group field cosmology.
Furthermore,  the supersymmetric loop quantum gravity have been studied only at the second quantized level \cite{super}-\cite{supera},
and  supersymmetry has never been studied at the level of third quantized loop quantum gravity or loop quantum cosmology.
  In fact, no third quantized
theory, including group field theory, has ever been supersymmetrized.
 Supersymmetry has been an important ingredient in the study
of M-theory \cite{mtheory}-\cite{mtheory1} and many
phenomenological models beyond the standard model are based
on it \cite{standard}-\cite{standard1}. It also provides
the an interesting
candidate for dark matter \cite{dark}-\cite{dark1}.
 Hence, it will be interesting to incorporate it in group
 field cosmology at a
third quantized level.  Motivated from this fact we also
construct a supersymmetric group field theory
with gauge symmetry in this paper.

\section{Loop Quantum Cosmology}
In loop quantum cosmology, the integration  is restricted to a fixed three dimensional cell of finite volume $ {\cal{V}}_0$.
 The integration is restricted in order to regulate the infinities coming from homogeneous fields
 being integrated over a non-compact space-time \cite{qazs}-\cite{qazr}.
The main  variables in   loop quantum gravity are  Ashtekar-Barbero connection $A^i_a $
and $E^a_i $.  In the loop quantum gravity the holonomies of the connection are defined as operators.
We will analyse the loop quantum cosmology for a massless scalar field $\phi$ in a spatially flat, homogeneous, and isotropic universe with
metric,
\begin{equation}
 ds^2 = -N^2(t) dt^2 + a^2 (t) \delta_{ab} dx^a dx^b.
\end{equation}
Here  $a(t)$ is the scalar factor and $N(t)$ is the lapse function.
Now $A^i_a = \gamma (\omega^i_0)_a,  $ where $\gamma$ is the Barbero-Immirzi parameter and $(\omega^i_0)_a$ is the Levi-Civita connection.
In loop quantum cosmology $A^i_a$  only depends on $c = \pm {\cal{V}} _0^{1/3} N^{-1} \gamma a'$,
 with $a'$ being the time derivative of $a$.
Furthermore, the  variable promoted to operators are $p = \pm a^2 {\cal{V}} _0^{2/3}$,  and $exp(i\mu c)$, where $c$ is conjugate to $p$, and
 $\mu$ is a function of $p$.
In the improved dynamics scheme \cite{a}, $\mu = |p|^{-1/2}$, and this forms a basis of the eigenstates of the volume operator
$\cal{V}$ with
\begin{equation}
 {\cal{V}} |\nu \rangle = 2 \pi \gamma G |\nu| |\nu \rangle,
\end{equation}
where $\nu = \pm a^2 {\cal{V}}_0 /2\pi \gamma G$ has the dimensions of length.
The variable conjugate to $\nu$ is $b$ and so
 the main variable for the gravitational sector  are
$\nu$ and $\exp( i \lambda b)$, where $\lambda$ is a constant.
For the matter sector one choses the usual Schroedinger quantization scheme and
the tensor product of the kinematic Hilbert space of the gravitational sector
with the kinematic Hilbert space of the matter sector give the total Hilbert space of the theory.

For the dynamics of the theory the curvature of $A^i_a$ is expressed through the holonomy around a loop. The area of such a loop cannot
cannot be smaller than a fixed minimum area because
the smallest  eigenvalue of the area operator in loop quantum gravity is nonzero.
In Planck units the Hamiltonian constrain for a homogeneous isotropic universe with a massless scalar field $\phi$,  can thus be written as
\begin{equation}
 K ^2\Phi(\nu, \phi) = [E^2 - \partial^2_\phi] \Phi(\nu, \phi) =0,
\end{equation}
where we have  defined $ E^2 $ to be
\begin{eqnarray}
 E^2 \Phi(\nu, \phi) &=& - [B(\nu)]^{-1}C^+(\nu) \Phi(\nu+4 , \phi)
 - [B(\nu)]^{-1}C^0(\nu)\Phi(\nu, \phi)  \nonumber \\ && - [B(\nu)]^{-1}C^-(\nu)\Phi(\nu-4, \phi).
\end{eqnarray}
Here we have set $\nu_0 =4$.
For the usual choice of gauge in  loop quantum cosmology, we have
\begin{equation}
 C^-(\nu) =C^+(\nu- 4 ),
\end{equation}
and so we need to only specify the functions  $C^+(\nu), C^0(\nu)$ and $B(\nu)$.
These functions depend of the choice of the lapse function and the quantization scheme choose,
 for example with $N =1$  in improved dynamic scheme, we have \cite{a}
\begin{eqnarray}
C^+(\nu)&=&\frac{1}{12\gamma\sqrt{2\sqrt{3}}}\left|\nu+2\right|\left| \left|\nu+1\right|-
\left|\nu+3\right| \right|,\nonumber
\\C^0(\nu)&=&-C^+(\nu)-C^+(\nu-4),\nonumber
\\
B(\nu)&=&\frac{3\sqrt{2}}{8\sqrt{\sqrt{3}}\pi\gamma G}|\nu| \left| \left|\nu+1\right|^{\frac{1}{3}}
-\left|\nu-1\right|^{\frac{1}{3}} \right|^3.
\end{eqnarray}
 Similarly, in solvable loop quantum cosmology with $N =a^3$, we have \cite{as}
\begin{eqnarray}
C^+(\nu)=\frac{\sqrt{3}}{8\gamma}\left(\nu+2\right),&& C^0(\nu)=-\frac{\sqrt{3}}{4\gamma}\nu,
\nonumber\\
B(\nu)=\frac{1}{\nu}. &&
\end{eqnarray}
There is  also a super-selection for  subspace of the  form $| \nu +4n\rangle$, for all
integers  $n$ and  some state $|\nu\rangle$.
It may also be noted that the matter variable appear as the time variable in  loop quantum cosmology.
\section{Group Field Cosmology}
Just like the wave function of first quantized theories is viewed as a classical field in second quantized formalism, the wave function of the
second quantized theory should be viewed as a classical field in third quantized formalism. Hence, the wave function of loop quantum cosmology
will now be viewed as the classical field of group field cosmology.
So, we want a free bosonic field theory such that its classical equation of motion reproduce the Hamiltonian
constraint for loop quantum gravity.
The classical action for this bosonic field theory can be written as  \cite{m}
\begin{equation}
 S_{b} = \sum_\nu \int d\phi \, \,    \Phi (\nu, \phi)  K ^2 \Phi(\nu, \phi).
\end{equation}
As $K^2$ is not diagonal in $\nu$, it is useful to expand $\Phi(\nu, \phi)$ as
\begin{eqnarray}
 \Phi (b, \phi) &=& \sum_\nu \Phi(\nu, \phi)  exp ( i \nu b), \nonumber \\
\Phi ( \nu, \phi) &=& \frac{2}{\pi} \int_0^{\pi/2} db  \,\,  \Phi(b, \phi) exp (-i\nu b).
\end{eqnarray}
Now the bosonic action is given by
\begin{equation}
 S_b = \frac{1}{\pi} \int d\phi \, \,  \int db  \,\,  \overline{\Phi}(b, \phi) \tilde K^2 \Phi(b, \phi),
\end{equation}
where $\Phi(v, \phi) = \overline{\Phi} (\pi/2 -b , \phi)$ and
the inverse of $ \tilde K^2$ is the propagator of the theory.
To construct the Fock space, we expand these fields into modes which we promote to operators \cite{m}.
Each of these will create or annihilate universes.
The vacuum state is defined as the state annihilated
by all annihilation  operator,
$
 a_k |0\rangle =0,
$
and the Fock space is constructed as by the action of creation operators on this vacuum state. The vacuum state here corresponds to a state with no
geometry, matter field and topological structure. The topology, geometry and matter is created by the action of
creation operators on this vacuum state. Interactions  will now correspond to an interaction of these universes.
Hence, they will represent the process of splitting of or joining of universes.  Thus the topology changing
processes can be studied by including interaction terms to this bosonic free field theory.
 Each interaction will also correspond to  a conservation of a different quantity during the process of interaction of these universes.
So, for example the interaction potential of the form $\lambda \Phi^3( \nu, \phi)/3$, will correspond to a conservation of
 both $b$ and $p_\phi$, because
it will implement locality in both the scalar factor and the scalar field.

There is no reason to assume that the statistics of the theory
are bosonic and so we
can also construct a fermionic model for group field cosmology.
The choice of fermionic statics will change
how distribution of different universes in the multiverse.
In fact, attempts to explain the cosmological constant
 from a bosonic distributions of universes in the multiverse
 have led to a wrong value for the cosmological constant
\cite{cos1}. It will be interesting to perform
 a similar analysis with a fermionic distribution of universes
and try to calculate the value of the cosmological constant.
 This is one of the motivations for studding a fermionic
distribution of universes in the multiverse.

Fermionic statics
will not change the dynamics of a single universe because the fermionic
theory  will also be required to generate
the Hamiltonian constraint of loop quantum gravity.
Thus, in analogy with a fermionic
field in two dimensions \cite{susy},
  we define a fermionic field in group field cosmology as
 $\Psi_a (\nu, \phi) = (\Psi_1(\nu, \phi), \Psi_2(\nu, \phi) )$.
The spinor indices are
 raised and lowered by the second-rank
antisymmetric tensors $C^{ab}$ and $C_{ab}$, respectively,
\begin{eqnarray}
 \Psi^a (\nu, \phi)&=&  C^{ab}\Psi_b(\nu, \phi), \nonumber \\
 \Psi_a(\nu, \phi) &=& \Psi^b (\nu, \phi)C_{ba}.
\end{eqnarray}
These second-rank
antisymmetric tensors also satisfy $C_{ab}C^{cb} = \delta^c_a$.
So, we have $(\gamma^\mu)_{ab} = (\gamma^\mu)_a^c C_{cb}
= (\gamma^\mu)_{ba}$
\cite{bsusy}.
Now we
define $K_\mu = (E, \partial_\phi) $ and $\eta^{\mu\nu} = (1, -1) $.
Here $K_\mu$ acts
like $\partial_\mu$ of a ordinary two dimensional field  theory,
so in analogy with   $\partial_{ab} = (\gamma^\mu \partial_\mu)_{ab}$
and $\partial^2 = \partial^{ab}\partial_{ab}/2$ \cite{susy},
we also define  $ K_{ab} = (\gamma^\mu )_{ab}K_\mu $ and
\begin{eqnarray}
 K^2 &=& \frac{1}{2} K^{ab} K_{ab} \nonumber \\
&=& \frac{1}{2} (\gamma^\mu )^{ab}K_\mu(\gamma^\nu )_{ab}K_\nu
 \nonumber \\
&=& \eta^{\mu \nu } K_\mu K_\mu
 \nonumber \\
&=& [E^2 - \partial^2_\phi],
\end{eqnarray}
where we have used $\{\gamma^\mu, \gamma^\nu\} = 2 \eta^{\mu\nu}$.
Now we can write  the action for the fermionic group field cosmology
 in analogy
with the regular fermonic action in two dimensional
spinor formalism
by using the correspondence between
$K_\mu$ and $\partial_\mu$ \cite{susy},
\begin{eqnarray}
 S_{f} &=& \sum_\nu \int d\phi \, \, C^{bc} \Psi_c (\nu, \phi)
 (\gamma^\mu )_{b}^a K_\mu \Psi_a(\nu, \phi)\nonumber \\
&=& \sum_\nu \int d\phi \, \,  \Psi^b(\nu, \phi) K_b^a \Psi_a(\nu, \phi).
\end{eqnarray}
The classical equations of motion obtained from this action are
\begin{equation}
 K^b_a \Psi_b (\nu, \phi) =0.
\end{equation}
Now action on it by $K^c_b$, and
using $\{\gamma^\mu, \gamma^\nu\} = 2 \eta^{\mu\nu}$, we  get
\begin{equation}
 K^2 \Psi_c (\nu, \phi) =0.
\end{equation}
We can also construct the Fock space for this fermionic theory by first  expand these fields into modes and then promoting those modes  to operators.
Just like the bosonic theory, the vacuum state of the fermionic theory will also be annihilated by all annihilation operators
$
 a_k |0\rangle =0,
$
and the Fock space will be constructed  by the action
of creation operators on this vacuum state.
 This vacuum state again corresponds to a state with no
geometry, matter field and  topological structure.
Here the again modes denoted by $k$ correspond
to the  momentum $p_\phi$, which is the momentum
associated with the scalar field $\phi$ \cite{m}.
We can also introduce interactions for this fermionic theory or
an combination of both the
fermionic and bosonic theories.
\section{Supersymmetry}
In this section we will derive a supersymmetric generalization of the group field cosmology.
We will in fact construct a  super-group field  theory, such that the kinetic parts of its  bosonic  and fermionic component fields will
generate the Hamiltonian constraint of loop quantum gravity.
 In order to
do that  we first introduce $\theta_a$ as
two component anti-commuting parameters with odd Grassmann parity and let
$
 \theta^2 =  \theta_a C^{ab} \theta_b/2 =  \theta^a \theta_a/2.
$
Now  using this Grassman variable $\theta_a$, we define the generators of ${\cal_{N}} =1$ supersymmetry as
\begin{equation}
 Q_a = \partial_a - K_a^b \theta_b.
\end{equation}
The covariant derivative that commutes with this generator of supersymmetry can be written as
\begin{equation}
 D_a = \partial_a + K^b_a \theta_b,
\end{equation}
where
\begin{equation}
 \partial_a = \frac{\partial}{\partial \theta^a }.
\end{equation}
Hence, the full supersymmetric algebra between $D_a$ and $Q_a$ is given by
\begin{eqnarray}
 \{Q_a, Q_b\} &=& (\partial_a - K_a^c \theta_c) (\partial_b - K_b^d \theta_d)  \nonumber \\
 &&  + (\partial_b - K_b^d \theta_d) (\partial_a - K_a^c \theta_c)\nonumber \\ &=&  2 K_{ab}, \nonumber \\
\{ D_a, D_b\} &=& (\partial_a + K_a^c \theta_c) (\partial_b + K_b^d \theta_d)  \nonumber \\
 &&  + (\partial_b + K_b^d \theta_d) (\partial_a + K_a^c \theta_c)\nonumber \\ &=& - 2 K_{ab}, \nonumber \\
\{Q_a, D_b\} &=& (\partial_a - K_a^c \theta_c) (\partial_b + K_b^d \theta_d)  \nonumber \\
 &&  + (\partial_b +  K_b^d \theta_d) (\partial_a - K_a^c \theta_c)\nonumber \\ &=&  0,
\end{eqnarray}
because $(\gamma^\mu K_\mu)_{ab} = (\gamma^\mu K_\mu)_a^c C_{cb} = (\gamma^\mu K_\mu)_{ba}$ \cite{bsusy}.
Now we will write a supersymmetric theory in $ {\mathcal{N}}= 1$ superspace formalism so that it has manifest  $ {\mathcal{N}}= 1$ supersymmetry.
To do that, we define a superfield $\Omega( \nu, \phi, \theta)$, such that
\begin{equation}
 \Omega( \nu, \phi, \theta) = \Phi( \nu, \phi) + \theta^a \Psi_a( \nu, \phi) - \theta^2 F ( \nu, \phi).
\end{equation}
 Now we can write
\begin{eqnarray}
\Phi( \nu, \phi) = [\Omega( \nu, \phi, \theta)]_|, && \Psi_a( \nu, \phi)=  [D_a \Omega( \nu, \phi, \theta)]_|, \nonumber \\
F ( \nu, \phi)= [D^2 \Omega( \nu, \phi, \theta)]_|,
\end{eqnarray}
where $'|'$ means we set $\theta_a =0$ at the end of calculation.
The supersymmetric transformations $ \delta_{s}$ generated by $Q_a$ can also be calculated in analogy with the regular supersymmetric theory,
\begin{eqnarray}
 \delta_{s} \Phi( \nu, \phi) &=& - \epsilon^a ( \nu, \phi) \Psi_a ( \nu, \phi), \nonumber \\
\delta_{s} \Psi_a ( \nu, \phi) &=&- \epsilon^b( \nu, \phi) [C_{ab} F ( \nu, \phi)  + K_{ab}\Phi( \nu, \phi) ], \nonumber \\
\delta_{s} F ( \nu, \phi)&=&-\epsilon^a ( \nu, \phi)  K_a^b \Psi_b ( \nu, \phi).
\end{eqnarray}
These transformations are similar to the regular $\mathcal{N} =1$ transformations in three dimensions \cite{susy},
with the spacetime replaced by $(\nu, \phi)$ and $\partial_{ab}$ replaced by $K_{ab}$.
Now we can write a superfield theory in
this  $\mathcal{N} =1$  superspace. As the field theory will be written in $\mathcal{N} =1$, it will have manifest supersymmetry.
However, we have two operators $E$ and
$\partial_\phi$, so we have to work in analogy with a two dimensional superfield theory. It is well known that
$\mathcal{N} =1$ supersymmetry in three dimensions corresponds to $\mathcal{N} =2$ supersymmetry in two dimensions. This is because
by using the projection  $P_{\pm } = (1 \pm \gamma^3)/2$,
we can split the super-charge $Q_a$  as
$\epsilon^a Q_a = \epsilon^- ( \nu, \phi) Q_- + \epsilon^- ( \nu, \phi) Q_+$.
In two dimensions these $Q_+$ and $Q_-$ act as two independent super-charges \cite{bsusy}-\cite{susy1}. So,
our theory actually has $\mathcal{N} =2$ supersymmetry.
Thus, again in analogy with regular free supersymmetric theory, we write the following action,
\begin{eqnarray}
 S_{s} &=& \frac{1}{2}\sum_\nu \int d\phi \, \,
 [D^2 [\Omega ( \nu, \phi, \theta) D^2 \Omega ( \nu, \phi, \theta)]]_|
\nonumber \\
&=& \frac{1}{2} \sum_\nu \int d\phi \, \,
[ D^2 \Omega( \nu, \phi, \theta) D^2 \Omega( \nu, \phi, \theta)
 \nonumber \\ &&
 \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+ D^a\Omega( \nu, \phi, \theta) D_a D^2 \Omega( \nu, \phi, \theta) \nonumber \\ &&
 \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+ \Omega ( \nu, \phi, \theta) (D^2)^2 \Omega( \nu, \phi, \theta)]]_|
\nonumber \\ &=& \frac{1}{2}
 \sum_\nu \int d\phi \, \,
  [ F^2( \nu, \phi) + \Psi^b(\nu, \phi)
 K_b^a \Psi_a(\nu, \phi) \nonumber \\ &&
 \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,
+  \Phi (\nu, \phi)  K ^2 \Phi(\nu, \phi) ].
\end{eqnarray}
This is the required supersymmetric model whose bosonic part and fermionic parts generated the bosonic and fermionic
models of loop quantum cosmology.
We can now add interactions to this supersymmetric theory,
\begin{eqnarray}
  S_{si} &=& \sum_\nu \int d\phi \, \,  [D^2  [ f( \Omega( \nu, \phi, \theta))]]_| \nonumber \\
&=& \sum_\nu \int d\phi \, \,  [ f'' ( \Omega( \nu, \phi, \theta)) (D^a \Omega( \nu, \phi, \theta))^2 \nonumber \\ &&
 \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+ f'(\Omega( \nu, \phi, \theta)) D^2 \Omega( \nu, \phi, \theta)]_|   \nonumber \\
&=&    \sum_\nu \int d\phi \, \,  [f''(\Phi)( \nu, \phi) \Psi^2 ( \nu, \phi) + f'(\Phi)( \nu, \phi) F( \nu, \phi) ],
\end{eqnarray}
where $f( \Omega( \nu, \phi, \theta))$ is a general interaction term containing polynomial terms in $\Omega( \nu, \phi, \theta)$
 and $\Psi^2 ( \nu, \phi)= \Psi^a( \nu, \phi) \Psi_a( \nu, \phi)/2 $.
For example, the  superspace interaction term given by
\begin{equation}
f( \Omega( \nu, \phi, \theta)) = \frac{\lambda}{6} \Omega^3( \nu, \phi, \theta),
\end{equation}
generates  the following component action,
\begin{eqnarray}
 S_t &=& \frac{1}{2}\sum_\nu \int d\phi \, \,
[ F^2( \nu, \phi) +  \Psi^b(\nu, \phi) K_b^a \Psi_a(\nu, \phi)
\nonumber \\ &&
 \,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,
+  \Phi (\nu, \phi)  K ^2 \Phi(\nu, \phi)
+ 2\lambda \Phi (\nu, \phi)\Psi^2(\nu, \phi)  \nonumber \\ &&
  \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,
+ \lambda \Phi^2(\nu, \phi) F(\nu, \phi) ],
\end{eqnarray}
where $S_t =  S_{s} + S_{si}$.
It may be noted that the exact nature of the group field cosmology will effect the way individual universes interact in the multiverse.
As it is not possible at present to view other universes, it seems difficult to see which of these models is close
 to being the correct low energy description to some full quantum theory of gravity that describes the full multiverse. One way to deal with,
this problem is to view these models as describing the quantum states of our universe at Planck level and then trying to analyse it to see if
it can lead to some effective phenomenological description, at energy scales that we can measure.
So, it is  interesting consider these as patches  of the our universe that interact with one another.
Such a possibility has been discussed and used to analyse cosmological perturbations \cite{pat}-\cite{pat2}.
This possibility  has also been discussed in the framework of bosonic group field cosmology \cite{m}.
\section{Gauge Symmetry}
In analogy with regular particle physics, we can also incorporate gauge  symmetry into the
 third quantized group field cosmology. In fact, in this section,
 we
will constraint  supersymmetric group field cosmology  with gauge symmetry. To do that we first consider complex superfields
$\Omega(\nu, \phi, \theta) $ and $ \Omega^{\dagger}  (\nu, \phi, \theta)$, which are suitably contracted with generators of a Lie algebra,
  $\Omega(\nu, \phi, \theta) = \Omega^A(\nu, \phi, \theta) T_A,$ and $
\Omega^{\dagger}(\nu, \phi, \theta) = \Omega^{\dagger A}(\nu, \phi, \theta)  T_A$. Here
 $T_A$ are Hermitian generators of a Lie algebra
$[T_A, T_B] = i f_{AB}^C T_C$.
Now we define
\begin{eqnarray}
 \Omega^A(\nu, \phi, \theta) &=& \Phi^A( \nu, \phi) + \theta^a \Psi^A_a( \nu, \phi) - \theta^2 F^A ( \nu, \phi), \nonumber \\
 \Omega^{A \dagger}(\nu, \phi, \theta) &=& \Phi^{A \dagger}( \nu, \phi) + \theta^a \Psi^{A \dagger}_a( \nu, \phi)
 - \theta^2 F^{A \dagger} ( \nu, \phi),\nonumber \\
\Lambda^A (\nu, \phi, \theta) &=& \lambda^A (\nu, \phi)  + \theta^a  \mu_a^A (\nu, \phi) - \theta^2 \nu^A ( \nu, \phi).
\end{eqnarray}
These superfields transform under infinitesimal gauge
 transformations as
\begin{eqnarray}
  \delta \Omega^A(\nu, \phi, \theta) &=&  if_{CB}^A\Lambda^C (\nu, \phi, \theta)\Omega^B(\nu, \phi, \theta) ,\nonumber\\
\delta \Omega^{A \dagger}(\nu, \phi, \theta)  &=& -i f^A_{CB}\Omega^{C\dagger}(\nu, \phi, \theta) \Lambda^B(\nu, \phi, \theta).
\end{eqnarray}
These transformations can now be written as
\begin{eqnarray}
 \delta \Phi^A( \nu, \phi) &=&  i f_{CB}^A\lambda^C (\nu, \phi) \Phi^B( \nu, \phi), \nonumber \\
 \delta \Phi^{A \dagger}( \nu, \phi) &=&  - i f_{CB}^A (\nu, \phi) \Phi^{C\dagger}( \nu, \phi)\lambda^{B\dagger}, \nonumber \\
\delta\Psi^{A }_a( \nu, \phi) &=& if_{CB}^A [\lambda^C (\nu, \phi) \Psi^B_a( \nu, \phi) +\mu_a^C (\nu, \phi) \Phi^B( \nu, \phi) ],  \nonumber \\
\delta\Psi^{A \dagger}_a( \nu, \phi) &=& - if_{CB}^A [  \Psi^{C \dagger}_a( \nu, \phi)\lambda^{B \dagger}(\nu, \phi)
 + \Phi^{C \dagger}( \nu, \phi) \mu_a^{B \dagger} (\nu, \phi)],\nonumber \\
\delta F^A ( \nu, \phi) &=&  i f_{CB}^A [\lambda^C (\nu, \phi) F^B ( \nu, \phi) + \nu^C ( \nu, \phi)\Phi^B( \nu, \phi)
 \nonumber \\  && \,\,\,\,\,\,\,\,\,\,\,\, + \mu^{C a} (\nu, \phi)\Psi^B_a( \nu, \phi)], \nonumber \\
\delta F^{A \dagger} ( \nu, \phi) &=& - i f_{CB}^A [ F^{C\dagger} ( \nu, \phi) \lambda^{B\dagger} (\nu, \phi)
 + \Phi^{C\dagger} ( \nu, \phi)\nu^{B\dagger} ( \nu, \phi)
 \nonumber \\  && \,\,\,\,\,\,\,\,\,\,\,\,+\Psi^{C a\dagger}( \nu, \phi) \mu_a^{B\dagger} (\nu, \phi)].
\end{eqnarray}
 Now the  super-derivative,  given by
$
 D_a = \partial_a + K^b_a \theta_b,
$
of these superfields does not transform
 like the original superfields. But
 we can define a super-covariant derivative for  these superfields
by requiring it   to transform like the original superfields.
 Thus, we obtain the following expression for the
 super-covariant derivative of these superfields
\begin{eqnarray}
  \nabla_a  \Omega^A(\nu, \phi, \theta)&=& D_a\Omega^A(\nu, \phi, \theta) -i f_{CB}^A\Gamma^C_a(\nu, \phi, \theta) \Omega^B(\nu, \phi, \theta),\nonumber\\
\nabla_a \Omega^{A \dagger}(\nu, \phi, \theta)  &=& D_a \Omega^{A \dagger}(\nu, \phi, \theta)
+ i  f_{CB}^A\Omega^{C \dagger}(\nu, \phi, \theta)  \Gamma^B_a (\nu, \phi, \theta).
\end{eqnarray}
If  the spinor superfield $\Gamma_a^A (\nu, \phi, \theta) $ is made to transform under
gauge transformations as
\begin{equation}
 \delta \Gamma^A_a(\nu, \phi, \theta)  = \nabla_a \Lambda^A(\nu, \phi, \theta),
\end{equation}
then the  super-covariant derivative of the
 scalar superfields $\Omega^A(\nu, \phi, \theta)$ and $\Omega^{A \dagger}(\nu, \phi, \theta)$ transforms
 under gauge   transformations like the original fields,
\begin{eqnarray}
  \delta \nabla_a \Omega^{A \dagger}(\nu, \phi, \theta)&=&  if_{CB}^A\Lambda^C (\nu, \phi, \theta) \nabla_a\Omega^{B \dagger}(\nu, \phi, \theta) ,\nonumber\\
\delta \nabla_a\Omega^{A \dagger}(\nu, \phi, \theta) &=& -i f_{CB}^A \nabla_a\Omega^{C \dagger}(\nu, \phi, \theta)\Lambda^B(\nu, \phi, \theta).
\end{eqnarray}

Now we can use these superfields to write an action for  supersymmetric group field cosmology with gauge symmetry.
To that we first define $\Gamma_a(\nu, \phi, \theta)$ to be a matrix valued spinor
 superfield which is suitable contracted with generators of a Lie algebra,  $ \Gamma_a(\nu, \phi, \theta) = \Gamma_a^A(\nu, \phi, \theta) T_A $.
Now define the components of this superfield  $\Gamma_a(\nu, \phi, \theta)$ to be
\begin{eqnarray}
A_a (\nu, \phi)&= &[\Gamma_a(\nu, \phi, \theta)]_|, \nonumber \\
 X(\nu, \phi)&= & -\frac{1}{2}[D^a \Gamma_a(\nu, \phi, \theta)]_|, \nonumber \\
A^\mu (\nu, \phi)&= &- \frac{1 }{2} [ D^a (\gamma^{\mu })_a^b \Gamma_b (\nu, \phi, \theta)]_|, \nonumber \\
 X_a(\nu, \phi) &= &\frac{1 }{2}[D^b D_a  \Gamma_b(\nu, \phi, \theta)]_|. \label{csf}
\end{eqnarray}
It may be noted that the  gauge field $A^\mu (\nu, \phi)$ in group field  cosmology is analogy to ordinary gauge field.
Thus, $\Gamma_a(\nu, \phi, \theta)$ is the supersymmetric generalization of the  gauge field in group field cosmology.
We also
 define a field strength for this supersymmetric gauge field in group field cosmology as,
\begin{equation}
 \omega_a (\nu, \phi, \theta) = \nabla^b \nabla_a \Gamma_b(\nu, \phi, \theta).
\end{equation}
Now we can write the action for the gauge theory as
\begin{eqnarray}
S_{ga} &= &\sum_\nu \int d\phi \, \,    [D^2 [\Omega^{\dagger} ( \nu, \phi, \theta) \nabla^2 \Omega ( \nu, \phi, \theta)
\nonumber \\ && \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + \omega^a ( \nu, \phi, \theta) \omega_a( \nu, \phi, \theta)  ]]_|.
\end{eqnarray}
As this theory has a gauge symmetry, we need to fix a gauge before we can quantize it. Thus, in analogy with
ordinary supersymmetric gauge theory \cite{susy}, we chose the following gauge,
\begin{equation}
 D^a \Gamma_a ( \nu, \phi, \theta) =0.
\end{equation}
This can be incorporated at a quantum level by adding the following gauge fixing term to original classical action,
\begin{equation}
 S_{gf} = \sum_\nu \int d\phi \, \,    [D^2[ B ( \nu, \phi, \theta) D^a \Gamma_a ( \nu, \phi, \theta)]]_|,
\end{equation}
where $B( \nu, \phi, \theta) = B^A( \nu, \phi, \theta) T_A$ is a
auxiliary superfield. This field can be integrated out, in the path
integral formalism, to impose the
gauge fixing condition. Now we can write the ghost term by  gauge transforming the gauge fixing condition, replacing
$\Lambda(\nu, \phi, \theta) = \Lambda^A (\nu, \phi, \theta) T_A$
 by ghost superfield $ {C}( \nu, \phi, \theta) =  {C}^A( \nu, \phi, \theta)  T_A$
and then  contracted it with anti-ghost superfield $ \overline{C}( \nu, \phi, \theta) = \overline{C}^A( \nu, \phi, \theta)  T_A$ \cite{fdfd},
\begin{equation}
  S_{gh} = \sum_\nu \int d\phi \, \,
 [D^2[\overline C( \nu, \phi, \theta) D^a \nabla_a C
( \nu, \phi, \theta)]]_|.
\end{equation}
Thus, the sum of the  ghost and gauge fixing terms is given by
\begin{eqnarray}
 S_{gh} + S_{gf} &=& \sum_\nu \int d\phi \, \,    [D^2[ B ( \nu, \phi, \theta) D^a \Gamma_a ( \nu, \phi, \theta) + \nonumber \\
&& \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,
\overline C( \nu, \phi, \theta) D^a \nabla_a C ( \nu, \phi, \theta)]]_|.
\end{eqnarray}
It may be noted that in analogy with ordinary supersymmetric gauge theory the a mixing of $D_a$ and $\nabla_a$
occurs in the ghost term.
 Now  we can write the vacuum to vacuum transition as
\begin{equation}
 Z = \int D M \,  exp (iS),
\end{equation}
where $DM =  D\Omega^{\dagger } D \Omega DC D\overline{C} DB D\Gamma_a $ and
$ S = S_{ga}+  S_{gh}+  S_{gf}$, without setting $\theta_a = 0$. Now we can derive the propagators for this theory and the Feynman
 rules using standard field theory techniques.
What is
interesting is that this action allows for  processes like  a bosonic universe with the wave function $\Phi( \nu, \phi) $ or a
fermionic universe with the wave function $\Psi(\nu, \phi)$ to be interact with each other via exchange
of third quantized  virtual gauge universe $A_\mu (\nu, \phi)$. Thus, the big bang in this model can be viewed as creation of a pair of universe
by the interaction of two other universes via the exchange of some third quantized virtual gauge universe.
Just like the first quantized wave function is
viewed as a classical field in second quantization, the second quantized wave function can be viewed as a classical
field in third quantization. Now, we know that the particle creation can only be consistently explained
in the framework of a second quantized theory. From the perspective of a first quantized theory the wave function of
a particle will suddenly disappear when it get annihilated. Similarly, the wave function of a particle will suddenly appear
when it forms. However, from the perspective of a second quantized theory when the first quantized wave function is
viewed as a classical field, the disappearance or the appearance of a particle can be easily understood in the light of an
interactive field theory. Similarly,  from the perspective of second quantization the wave function of the universe suddenly appears at the big bang.
Similarly, the wave function of a black holes seems to disappear during the process of evaporation of a black hole. However,
from a third quantized perspective, it is just the creation and annihilation of these spacetimes and matter configurations, that can
be understood in the light of interactive group field cosmology.

In order to consistently analyse the supersymmetric gauge field cosmology we need to analyse the third quantized
BRST symmetry of the theory. These symmetries can be use to show that various topology changing
processes are unitarity from a third quantized perspective.  So,  we will analyse the third quantized BRST  symmetry for this theory.
The sum of the ghost term with the  gauge fixing term can  be written as a total third quantized BRST variation,
\begin{equation}
  S_{gh} + S_{gf} = \sum_\nu \int d\phi \, \,    s \, [D^2 [\overline{C}  ( \nu, \phi, \theta) D^a \Gamma_a ( \nu, \phi, \theta) ]]_|,
\end{equation}
where the third quantized BRST transformations are given by
\begin{eqnarray}
  s\,\Omega^A(\nu, \phi, \theta)&= &if_{CB}
  ^AC^C (\nu, \phi, \theta)\Omega^B(\nu, \phi, \theta) , \nonumber \\
s\,  \Omega^{A \dagger}(\nu, \phi, \theta)  &= & -i f^A_{CB}\Omega^{C\dagger}(\nu, \phi, \theta) C^B(\nu, \phi, \theta),
\nonumber \\
s\, C^A( \nu, \phi, \theta)&= &f^A_{CB}C^{C\dagger}(\nu, \phi, \theta) C^B(\nu, \phi, \theta),
\nonumber \\
 s\, \Gamma^A_a ( \nu, \phi, \theta) &= & \nabla_a C^A( \nu, \phi, \theta),
\nonumber \\
s\, \overline{C} ^A( \nu, \phi, \theta) &= & B^A( \nu, \phi, \theta),
 \nonumber \\
 s\, B^A( \nu, \phi, \theta) &= &0.
\end{eqnarray}
Here we have  obtained third quantized BRST by following the
procedure for obtaining the   BRST transformations of any gauge field theory \cite{BRST}-\cite{BRST1}.
 Thus, the third quantized BRST transformations of the gauge and matter fields is there gauge transformations
with the gauge parameter replaced by the ghost superfield. The third quantized BRST transformation of the
ghost superfield is given by by contacting two ghost
superfields with the structure constant of the gauge group and the third quantized BRST transformation of anti-ghosts is
the auxiliary superfield used to impose
gauge fixing condition. Finally, The third quantized BRST transformation of the auxiliary field vanishes.
 These third quantized BRST
transformations are nilpotent, $s^2 =0$, and so the sum of the ghost term with the gauge fixing term is invariant under these third quantized
BRST transformations.
As the third quantized BRST transformations of the original third quantized classical action is only a ghost valued gauge transformations,
it is also invariant
under these third quantized BRST transformations. Thus, the full action is invariant under these transformations, $
 s\, S  =0$. These third quantized BRST transformation can be used to construct a third quantized BRST charge and analyse
 the unitarity of the supersymmetric group field
cosmology with gauge symmetry.

\section{Conclusion}
In this paper we have analysed the group field cosmology.
 To do so, we have analysed the loop quantum cosmology in third quantized formalism. Thus, after reviewing the
basics of loop quantum cosmology, we constructed both
 bosonic and fermionic field theories   whose classical equations of motion  generated the
the Hamiltonian constraint of loop quantum cosmology. We  discussed the effect of interactions for these models of group field cosmology.
We also constructed a supersymmetric theory in $ {\mathcal{N}} = 1 $ superspace
formalism with gauge symmetry.
 The component fields of this supersymmetric gauge theory  satisfied the Hamiltonian constraint of loop quantum
cosmology.
We  analysed the interactions in this supersymmetric  theory and derived the corresponding interactions for the bosonic and
fermionic fields component fields theories. In the supersymmetric gauge theory, two real universes interacted by an exchange of virtual universes.
The big bang in this model was  viewed as creation of a pair of universe
by the interaction of two other universes via the exchange of some third quantized virtual gauge universe.
We also analysed the third quantized BRST of this supersymmetric group field theory with gauge symmetry.
 We stress the fact that we did not
supersymmetrize the theory at the level of second quantized fields but at
the level of third quantized fields.
This is the first time that supersymmetry has been studied at the level of a third quantized field theory.

It would be interesting to carry to analyse further. One direction to carry out this analyses would be to perform it for
the full group field theory. Both the fermionic and bosonic group field theory have been studied and color has also been
incorporated in group field theory \cite{color}-\cite{colora}. It would thus be interesting to
analyse a supersymmetric group field theory with gauge symmetry.
 We will be able to develop that theory in analogy with the usual gauge theory.
It will be possible to derive this  theory in superspace formalism with manifest supersymmetry.
 This analysis can also be initially done in the ${\mathcal{N}}= 1$
superspace formalism. however, it will  be nice to analyse this present theory or the full group theory in the
extended superspace formalism, with higher amount of manifest supersymmetry.
It would also be interesting to analyse string field theory via group field cosmology.
 String theory can be viewed as two dimensional gravity
coupled to matter fields. Thus, it can be quantized via loop
quantum gravity. After that a group field cosmology of string theory
can be constructed. This
way we will be able to incorporate higher order string interactions in
string field theory. In fact, the   Hamiltonian constraint for
string in loop quantum gravity has already been analysed
\cite{stringfi}.
Now this
 Hamiltonian constraint string in loop quantum gravity,
$K^2_{string}\Phi =0$,  can view this as a kinetic part
of a group field theory action, $\Phi K^2_{string} \Phi$.
Then interactions can and introduce  through polynomial
interaction $\Phi^n$.
Hence, we will be able to go beyond the regular
cubic string field theory by using group field cosmology \cite{stringfi1}.
In fact, we can use all the results of this paper to analyse
this string field theory.
We could also try to construct the a superstring field theory of
 heterotic string as a group field
theory with the symmetry group $E_8\times E_8$ or $Spin(32)/Z_2$.
We can then incorporate arbitrary interactions for these string field
theories
via the standard methods of group field theory.


\section{Introduction}
The nuclear electric polarizability $\alpha_E$ is related to the response of a nucleus to an externally applied electric field. It is an interesting observable because it encapsulates information
about the excitation spectrum of a nucleus. Recently, it has attracted a lot of
attention both for light nuclei, see {\it e.g.}~\cite{Moro}, and for heavy
nuclei, see {\it e.g.}~\cite{Pb_RCNP}.
For light systems  the nuclear polarizability is relevant
in the extraction of nuclear quantities from atomic spectroscopic
measurements. The atomic energy levels are affected by  polarization of the
nucleus due to the electric field of the surrounding electrons. Such nuclear
structure correction, which is proportional to
$Z^3\alpha_E/a_0$ \cite{noert} where $a_0$ is the Bohr radius, needs to be
considered in the sophisticated quantum electrodynamics  calculations of the
atomic levels that  allow the extraction of charge radii from isotope shifts
measurements of unstable nuclei, see Refs.~\cite{Wan04, Maxime} and
\cite{Mue07} for $^6$He and $^8$He, respectively.
An even larger effect of the nuclear structure correction coming from the polarizability is expected in muonic atoms, as the  muon mass
is larger than the electron mass and the orbital radius is smaller.
This will be relevant to the proposed $\mu ^4$He and  $\mu ^3$He experiments \cite{rolf_exp} that aim at measuring the nuclear charge radius of $^4$He and $^3$He from the Lamb-shift, to be compared to electron scattering data.

 The nuclear electric polarizability of Helium isotopes is interesting for the several above mentioned reasons.
It has been already directly measured or extracted from experimental data for
the $^{3,4,6}$He isotopes \cite{Moro}. In the case of $^3$He it is worth
mentioning that the data from elastic scattering on Pb at energies below the
Coulomb barrier  \cite{Goe91} are in disagreement with estimates based on
calculations of the photo-absorption cross section~\cite{Moro}, the latter
being about a factor of 2 smaller. It is also worth noticing that the
theoretical calculations are in agreement with photo-absorption
experiments, and that the band
spanned by using different Hamiltonians in the calculations is smaller than
the difference between the data taken from photo-absorption cross section and from
ion scattering experiments. The data analysis involved in the  latter approach
is quite delicate, because one has to separate effects of the nuclear force
from Coulomb effects.

In Ref.~\cite{Stetcu} the polarizabilities of several Hydrogen and Helium isotopes were calculated with  {\it ab-initio} methods.
 Among the Helium isotopes, $^3$He and $^4$He were dealt with, but no prediction for $^6$He was provided. It is the aim of this paper to fill this gap.

The $^6$He is known as a  halo nucleus, where a tightly bound  $^4$He core is surrounded
by two neutrons~\cite{Tan96}.  It happens to be the lightest of the known halo nuclei
 and it is a Borromean nucleus, because the
two-neutron and neutron-core subsystems are unbound, but the
three-body system is held together.
Due to the very small separation energy which characterizes halo nuclei, one expects the
polarizability of $^6$He to be much larger than that of the tightly bound $^4$He isotope.
Experimental data indicate this behavior.
In this paper we would like to see whether microscopic calculations
reproduce the experimental values and lead to a result where
 $\alpha_E(^6{\rm He}) \gg \alpha(^4{\rm He}$).

We perform a microscopic study of the nuclear polarizability $\alpha_E$ for
$^6$He and compare it to $^4$He. We limit our study  to simple semi-realistic
two-body forces.
For that we use the  hyperspherical harmonics method with an effective interaction,
EIHH, to speed up the convergence~\cite{Nir,EIHH}.
 The polarizability is calculated as a sum rule of the dipole response function using the Lanczos algorithm and also integrating the photo-absorption cross section calculated with the Lorentz integral transform method~\cite{ELO94}.

This paper is organized as follows. In Section II we describe in details the theoretical calculation of the polarizability.
In Section III we present our results and in Section IV we make a comparison
with experiment. Finally, we conclude in Section V.

\section{Theoretical aspects}

The nuclear electric polarizability in the unretarded dipole approximation
is defined by the expression
\begin{equation}
\alpha_E=2 \alpha \sum_{n\ne0} \frac{|\langle n| {D_z}| 0\rangle|^2}{E_n -E_0}\,
\label{pol}
\end{equation}
where $\alpha$ is the fine structure constant, $D_z$ is the unretarded dipole
operator and  $E_{0/n}$ are the energies of the nuclear ground and excited states $|0\rangle$ and $|n\rangle$, respectively.
This observable is clearly related to the photo-absorption cross section and to the dipole
response function.
The photo-absorption cross section $\sigma_\gamma(\omega)$ of a nucleus
is given by
\begin{equation}
\sigma_\gamma(\omega)= 4\pi^2 \alpha
\, \omega\, R(\omega)\,,\label{sigma}
\end{equation}
where $R(\omega)$ is the response function.
In the unretarded dipole approximation
\begin{equation}
R(\omega) = \sum_{n,\bar{0}} |\langle n| D_z| 0\rangle|^2\,
            \delta(\omega-E_n+E_0) \,,
\end{equation}
where $\bar{0}$ indicates an average on the initial angular
momentum projections.
The dipole operator is given by $D_z=\sum_{i=1}^A z_i \tau_i^{3}/2$,
where $A$ is the number of nucleons and $\tau^3_i$ and $z_i$ are the third
component of the isospin
operator and the coordinate of the $i$th particle in the center of mass frame, respectively.
One can recover the expression for $\alpha_E$ in Eq.~(\ref{pol}) by calculating
sum rules of the photo-nuclear cross section.
The various  moments of $\sigma_{\gamma}$ are defined as
\begin{equation}
m_n(\bar\omega) \,\,\, \equiv \int_{\omega_{th}}^{\bar\omega}\,d\omega \,\omega^n\,
\sigma_\gamma(\omega)\,,\label{moments}
\end{equation}
where $\omega$ is the photon energy and $\omega_{th}$ and $\bar\omega$ indicate threshold energy and upper integration limit, respectively.
Assuming that $\sigma_\gamma(\omega)$ converges to zero
and utilizing the closure  of the  eigenstates of
the nuclear Hamiltonian $H$, one can relate the polarizability to the
  $n=-2$ sum rule,
\begin{eqnarray}
\alpha_E&=& \frac{m_{-2}(\infty)}{2\pi^2}
 = 2 \alpha
\,\sum_n \frac{|\langle n| D_z|0\rangle|^2}{E_n-E_0}\,.\label{mPSR}
\label{alphaEPSR}
\end{eqnarray}
The polarizability $\alpha_E$ can be calculated
with the Lanczos algorithm using a proper pivot. It is useful to rewrite Eq.~(\ref{alphaEPSR})
as
\begin{eqnarray}
\nonumber
\alpha_E&=& 2\alpha \langle 0| {D}^{\dagger}_z \frac{1}{H-E_0} D_z |0\rangle\,\\
\label{DHm1D}
&=& - 2\alpha \langle 0|D_z^{\dagger}D_z|0\rangle \langle \phi_0| \frac{1}{E_0-H} |\phi_0\rangle\,,
\end{eqnarray}
with
\begin{equation} \label{start}
  |\phi_0\rangle =\frac{D_z|0\rangle }{\sqrt{\langle
  0| D_z D_z|0 \rangle }} \,.
\end{equation}
Starting from the ``pivot'' of Eq.~(\ref{start})
where the ground state $|0\rangle$ is obtained by solving the Schr\"{o}dinger equation,
 $\alpha_E$ can be expressed as a continued fraction containing the Lanczos coefficients~\cite{Dagotto}
\begin{equation}
   a_i=\langle \phi_i |H|\phi_i \rangle \mbox{,}\,\,\,\,\,\,
b_i=\langle \phi_{i+1}|H|\phi_i \rangle \mbox{,}
\end{equation}
where the $|\phi_i\rangle$ form the Lanczos orthonormal basis $\{|\phi_i
\rangle \mbox{,}i=0,\ldots \}$. In fact one has
\begin{equation}
\alpha_E= - 2\alpha \langle 0|D_z^{\dagger}D_z|0\rangle \cfrac{1}{E_0-a_0- \cfrac{b_1^2}{E_0-a_1 -\cfrac{b^2_2}{E_0-a_3 \dots}}} \,.
\label{cfrac}
\end{equation}

In this work we calculate the polarizability in two different ways. On the one hand we
utilize Eq.~(\ref{cfrac}).
On the other hand we obtain $m_{-2}$ by integrating our results for the
 total photo-absorption cross section  calculated with
the Lorentz integral transform (LIT) method \cite{ELO94}.
In Ref.~\cite{He6Sonia1, He6Sonia2}
we have presented microscopic
calculations of the $^6$He $\sigma_{\gamma}$ with semi-realistic potential
models. Here we  use lager model spaces  which are nowadays available.
The LIT, an integral transform with a Lorentzian kernel, is defined as
\begin{equation}
{\cal L}(\sigma_R,\sigma_I) = \int d\omega {\frac {R(\omega)} {(\omega-\sigma_R)^2
                              + \sigma_I^2} } \,.
\end{equation}
The LIT is also typically calculated   using the Lanczos technique explained above, (see~\cite{Mario}). In fact it can be re-expressed as
\begin{equation} \label{loreins}
{\cal  L}(\sigma_R, \sigma_I)=-\frac{1}{\sigma_I}\,
            \langle 0 | D_z D_z|0\rangle\,
            \mbox{Im} \,\{\langle \phi_0 |
                      \frac{1}{z-H}|\phi_0\rangle\} \mbox{}
\end{equation}
with $z=E_0+\sigma_R+ i\sigma_I$.
It is evident that the LIT in (\ref{loreins})
is also a
 continued fraction as in Eq.~(\ref{cfrac}), where $E_0$ is replaced by a complex  $z=E_0+\sigma_R+ i\sigma_I$.
Once  ${\cal L}(\sigma_R,\sigma_I)$ is calculated, one can invert the LIT \cite{inversion}  to get
$R(\omega)$, and thus $m_{-2}$.
The two methods have to agree within the numerical uncertainty. However, with
the first method one avoids the complications
introduced by the inversion procedure.

Given the Hamiltonian $H$, the calculation of $\alpha_E$ in both ways is based on the EIHH~\cite{EIHH}
expansion of the wave function.
This  approach  is translational invariant, being constructed with the Jacobi coordinates.
We use different semi-realistic potential models for our calculations. Following
Ref.~\cite{He6Sonia1}, we will use the
 Minnesota (MN) potential \cite{Minnesota}
\begin{eqnarray}
\nonumber
V_{ij}&=&\left[V_R+\frac{1}{2}\left(1+P^{\sigma}_{ij}\right)V_T+\frac{1}{2}\left(1-P^{\sigma}_{ij}\right)V_S\right] \times \\
\label{mnu}
&&\left[\frac{1}{2}u+\frac{1}{2}\left(2-u\right)P^r_{ij}\right] \,,
\end{eqnarray}
where $P_{ij}^{\sigma,r}$ are spin and space-exchange
operators, $V_{R,T,S}$ are parameterized as linear combination of Gaussians of
the two-body relative distance and $u$ is a parameter.
This force reproduces the $S-$wave nucleon-nucleon phase shifts and correctly
binds the deuteron. It renormalizes effects of the tensor force into its
central component.
A typical value for $u$ in the Minnesota potential is $u=1$, as we used in~\cite{He6Sonia1}. Here we will explore the variation of this parameter by choosing $u\ge1$.
The mixing parameter $u$ does not affect the dominant $^1S_0$ and $^3S_1$
waves in the NN interaction, but only affects the $s=1$, $t=1$ channels, where
the dominant components are the $P-$waves ($^1P_1$ and $^3P_{0,1,2}$).
For $u=1$ there are no $P-$waves, they contribute only for $u>1$.
Thus, changing $u$ mostly affects $^6$He, without substantially changing $^4$He.
 Because in~\cite{He6Sonia1, He6Sonia2} we also used the Malfliet-Tjon (MTI-III) \cite{MTI-III},
 and the Argonne AV4' \cite{AV4'} potentials, we will present some results  with these interactions as well.
The Minnesota potential has been recently used in  a microscopic cluster model  calculation of $^6$He  \cite{Brida}  and in the Gamow shell model approach \cite{Papadim}  for
$^6$He and $^8$He.

\section{Results and Discussion}

\begin{figure}
\includegraphics[scale=0.35,clip=]{He4_polarizability_new.eps}
\caption{(Color online) $^4$He polarizability:  calculations with the Minnesota potential for two different $u$ values are shown as a function of the grand-angular momentum quantum number $K_{\rm max}$.
The polarizabilities obtained from realistic two- and three-body interactions \cite{Stetcu, Gazit} are  presented as a light (blue) band. Experimental data from \cite{Moro, Friar} are given by the dark (green) band.
}
\label{fig_He4}
\end{figure}

The main focus of this work is to study the $^6$He polarizability.
We start, however, the discussion with the $^4$He nucleus.
In Fig.~\ref{fig_He4}, we show the results of $\alpha_E$
calculated via the Lanczos coefficients, as in Eq.~(\ref{cfrac}).
The ground state $|0\rangle$ and the Lanczos ``pivot'' $ |\phi_0 \rangle$
are given in terms of the EIHH expansions.  While for the ground state the expansion is
characterized by an even hyperspherical grand-angular quantum number $K_{\rm max}$ and total
isospin $T=0$, $T_z=0$, $ D_z|0 \rangle$ has to be expanded on odd grand-angular quantum number $K'_{\rm max}$, where the isospin in the final state is $T'=1$.  Fig.~\ref{fig_He4} shows the convergence of $\alpha_E$ as a function of $K_{\rm max}$, where for each point $K_{\rm max}+1$ is used for the Lanczos ``pivot''.  We show our results for the Minnesota potential with $u=$1 and 1.20.
The convergence is very good, the dependence on $u$ is mild and the results are
very close to calculations where realistic nucleon-nucleon (NN) and three-nucleon (3N) forces have been used. For the latter one, results for effective field theory potentials were presented in \cite{Stetcu} leading to $\alpha_E=0.0683(8)(14)$ fm$^3$ (corresponding to the upper light (blue) band in Fig.~\ref{fig_He4}).
The error bar of this calculation is accounting for the convergence error $0.0008$ fm$^3$
 and also for the uncertainty in the underlying dynamics $0.0014$ fm$^3$.
We also show the  results  of   $\alpha_E$ for the Argonne $v_{18}$ two-body force and Urbana IX three-body force of Ref.~\cite{Gazit} leading to $\alpha_E=0.0655(4)$
 (corresponding to the lower light (blue) band in Fig.~\ref{fig_He4}), where the error bar comes from convergence only.
The experimental data are shown as a darker (green) band. These include the more recent evaluation of Ref.~\cite{Moro} based on the Arkatov {\it et al.}~\cite{Arkatov} experimental measurement of the photo-absorption cross section and an older result reported in Ref.~\cite{Friar}, based on earlier measurements by Arkatov {\it et al.}~\cite{Arkatov_old}. We would like to point out that the semi-realistic Minnesota potentials lead to a value of the polarizability which is consistent with realistic calculations and is only about 15$\%$ smaller than the average value in the experimental band.

We can also calculate the polarizability by integrating the photoabsorption cross section obtained with the LIT method. We get perfect agreement as with the Lanczos coefficients. For example, for the standard Minnesota potential where $u=1$ and for a $K_{\rm max}=12/13$ model space, the Lanczos method gives
$\alpha_E=0.06360$ fm$^3$ and integrating $\sigma_{\gamma}$ up to 120 MeV we get 0.06336 fm$^3$, with just a $0.4\%$ difference.

\begin{figure}
\includegraphics[scale=0.35,clip=]{Polarizability_vs_K.eps}
\caption{(Color online) The polarizability of $^6$He
as a function of the grand-angular momentum $K_{\rm max}$ for different semi-realistic interactions: the Minnesota potential with $u=1-1.2$ and the MTI-III potential.}
\label{fig_alphaE}
\end{figure}

Now we move to the $^6$He nucleus.
We first  calculate $\alpha_E$ from the Lanczos coefficients.
Also in this case the ground state is  expanded on even hyperspherical grand-angular quantum number $K_{\rm max}$, but the total
isospin  is $T=1$, $T_z=-1$, and $ D_z|0 \rangle$  is expanded on odd $K_{\rm max}'=K_{\rm max}+1$. In this case though, the final isospin can be
$T'=1$ or $T'=2$. This leads to two possible isospin channels that are calculated
separately and that open up at different energies. Experimentally, the $T=1$ channel opens up at photon energy $\omega_{th}=0.975$ MeV, while the $T=2$ channel opens up at $\omega_{th}=22.77$ MeV,  with $\gamma~^6$He $\rightarrow$ $^3$H$\,n\,n\,p$.
 Due to the inverse energy weight in  Eq.~(\ref{mPSR}),
the $T=2$ channel is expected to be less relevant to $\alpha_E$.
From our calculations we find that the percentage contribution of the T=2 isospin channel to the total polarizability changes from $2\%$ to $4\%$ when varying $u$ from 1 to 1.20 in the Minnesota potential.

\begin{table}
\caption{Results of the EIHH calculation with $K_{max}=12$ for the different $u$ value of the Minnesota potential. The values for the energies are in MeV.}
\label{Table1}
\begin{center}
\begin{tabular}
{c | lll}
\hline
Potential & $E_0$($^4$He) & $E_0$($^6$He) & $S_{2n}$($^6$He)\\
\hline\hline
Minnesota &  &  & \\
u=1.00 & -29.949 & -30.45 & 0.50\\
u=1.05 & -29.978 & -31.13 & 1.15\\
u=1.10 & -30.007 & -31.88 & 1.87\\
u=1.15 & -30.037 & -32.72 & 2.68\\
u=1.20 & -30.069 & -33.65 & 3.59\\
MTI-III & -30.760 & -32.24  & 1.48\\
\hline\hline
\end{tabular}
\end{center}
\end{table}

In Fig.~\ref{fig_alphaE}, we present a similar plot as in Fig.~\ref{fig_He4} for $^6$He with semi-realistic interactions.
 We observe a much slower convergence of $\alpha_E$ for $^6$He than for $^4$He with all the potentials employed.
By looking at the different $u$ values in the Minnesota potential, we see that
the convergence rate and the value of $\alpha_E$  substantially change  with $u$.
This is related to the variation of the binding energy, and consequently of the two-neutron separation energy, whose  numerical values are shown in Table~\ref{Table1} for completeness.
By increasing $u$ we are adding more $P-$wave interactions, which bring
additional binding to the $^6$He nucleus, while leaving $^4$He almost
unaffected. Naively, this makes $^6$He more tightly bound and thus more
difficult to polarize,  i.e., $\alpha_E$ gets smaller. For the value of $u=1$
the convergence of the polarizability is particularly slow, due to the fact
that $^6$He is barely bound, with  $S_{2n}=0.56$ MeV which is about a factor
of 2 smaller than the experimental value.
With the MTI-III potential we get a convergence pattern which is close to the Minnesota potential for $u=1.05$, because the prediction of $S_{2n}$ is similar with these two potentials, see Table~\ref{Table1}. The information that one gains from Fig.~\ref{fig_alphaE} is that by increasing $S_{2n}$ we can change the overall slope of the convergence pattern of $\alpha_E$.

\begin{figure}
\includegraphics[scale=0.35,clip=]{S2n_mncu.eps}
\caption{(Color online) $^6$He two-neutron separation energy as a
function of the grand-angular momentum $K_{\rm max}$ for the Minnesota potential and different $u$ parameters. The experimental value
is also shown.
}
\label{fig_S2n}
\end{figure}
For any considered value of $u$ though, it is clear that
our calculations reproduce the fact that the polarizability of the halo nucleus of $^6$He is much larger than that of the tightly bound $^4$He, the ratio being almost an order of magnitude.

Here we would like to point out that Brida and Nunes
in \cite{Brida}  have used the Minnesota potential with $u=1.15$ in a microscopic cluster model
and obtained a separation energy $S_{2n}=0.90(5)$ MeV.
This result is  different from the value we obtain and report
in Table~\ref{Table1}. Their calculation is performed without the Coulomb force, but its effect cancels in the separation energy. Because for $^4$He the value reported in \cite{Brida} for the binding energy is -30.85 MeV, which is in agreement with our value of -30.86(1) (with no Coulomb force), we think that the difference is due to the cluster assumption  made for $^6$He. We do not make such assumption and in convergence, the EIHH result is exact.
In  Fig.~\ref{fig_S2n}, we show that
the separation energy $S_{2n}$ is very well converged within the model space available for all
these potentials.

We can also calculate the polarizability by integrating the photoabsorption cross section obtained with the LIT method and compare it to the above results. We quote numbers for the  Minnesota potential with $u=1.05$ in the largest available model space $K_{\rm max}=12/13$.
 The Lanczos method gives
$\alpha_E=0.7542$ fm$^3$ and integrating $\sigma_{\gamma}$ up to 40 MeV (60 MeV) we get 0.7711 fm$^3$ (0.7827) fm$^3$.
Integrating the cross section we have a 3-4$\%$ difference, which is due
  to the fact that LIT is not completely converged and the inversion procedure
  introduces some numerical error.

\begin{figure}
\includegraphics[scale=0.35,clip=]{Polarizability_vs_S2n.eps}
\caption{(Color online) Correlation between $\alpha_E$ and
 $S_{2n}$ in $^6$He obtained with
the Minnesota potential varying the parameter $u$. The MTI-III and AV4' results
are also shown.
}
\label{fig_correlation}
\end{figure}

Now we would like to investigate the
dependence of the polarizability on the two-neutron separation energy. This can be achieved for example
 by plotting $\alpha_E$ versus $S_{2n}$ for the different values of
the parameter $u$ in the  Minnesota potential.
In Fig.~\ref{fig_correlation}, we can see that we find a correlation between $\alpha_E$ and  $S_{2n}$.
Calculations
have been performed with $K_{\rm max}=12$ ($K'_{\rm max}=13$). As an estimate of the theoretical error bar in the few-body method
we take  the difference between the largest possible calculation  with $K_{\rm max}=12$ and the
$K_{\rm max}=8$ result. We also present the data for the MTI-III and AV4'
potential (as used in \cite{He6Sonia2}) for completeness. The error bars for
the polarizability increase as the separation energy gets smaller. This is a
reflection of the slower convergence observed in Fig.~\ref{fig_alphaE}.
For the Minnesota potential $S_{2n}$ has a negligible error, hardly visible
in Fig.~\ref{fig_correlation}.
For the MTI-III and AV4' potentials the error in $S_{2n}$ is large because these interaction models are not as soft as the Minnesota force.

In Ref.~\cite{Stetcu} it was argued that the polarizability should roughly scale like the inverse square of the binding energy of a nucleus. For a halo system, like $^6$He, the relevant scale parameter
is the separation energy, rather than the binding energy. The $\alpha_E-S_{2n}$ dependence empirically observed in Fig.~\ref{fig_correlation} is compatible with such a behavior.

In order to reproduce the polarizability of a halo nucleus, it is expected that the halo structure, thus $S_{2n}$, should be correctly modeled, even if the absolute binding of $^4$He and $^6$He are not reproduced. Thus, one can estimate the value of $\alpha_E$ by choosing  $S_{2n}$ to be around the experimental value  and then calculate the corresponding polarizability. A value of $u$ that gives $S_{2n}$ close to experiment is $u=1.05$, where the convergence of $\alpha_E$ is slower than for larger values of $u$.
From a closer look to Fig.~\ref{fig_alphaE} and Fig.~\ref{fig_S2n} we can see that  also for $u=1.20$
the polarizability $\alpha_E$ is still increasing when $K_{max}$ becomes larger, even though
the separation energy is converged.
This means that the convergence of the polarizability is not only influenced by $S_{2n}$.
Other observable that is naturally related to the polarizability in the unretarded dipole approximation is the radius operator.

\begin{figure}
\includegraphics[scale=0.35,clip=]{alpha_vs_rskin.eps}
\caption{(Color online) Correlation between the nuclear electric polarizability of $^6$He and the skin radius with the Minnesota potential with different $u$. The four points for each $u$ value correspond to calculations with $K_{\rm max}=6,8,10$ and 12, moving from the left to the right.
}
\label{fig_alpha_skin}
\end{figure}

\begin{figure}
\includegraphics[scale=0.35,clip=]{r_skin_fit.eps}
\caption{(Color online) Neutron skin radius $r_{\rm skin}$  of $^6$He with the Minnesota potential and $u=1.05$, as a function of the grand-angular momentum quantum number $K_{\rm max}$.
The curve is a fit to the calculated points, used to extrapolate to infinite model space.
}
\label{rskin_fit}
\end{figure}

In  recent papers \cite{Pb_theory, Shlomo} the correlation between the
polarizability and the neutron skin of the $^{208}$Pb nucleus was studied
within the nuclear density functional theory framework.
In the following, we will  investigate the same correlation for $^6$He, even though $^6$He is a  different system. For halo nuclei, one refers to the halo radius, rather than the skin radius, but clearly the observable
\begin{equation}
r_{\rm skin}=r_n -r_p \,,
\end{equation}
where $r_n$ and $r_p$ are the mean point-neutron and point-proton radii, can be uniquely defined.
In our recent work~\cite{sonia_radii}, such observables have been calculated  for $^6$He from realistic two-body potentials in the EIHH method. Here, instead, we use the same semi-realistic interaction as for the $\alpha_E$ calculations.
 In Fig.~\ref{fig_alpha_skin}, we show a plot of $\alpha_E$ versus $r_{\rm skin}$
for different model spaces and for three different values of $u$ in the Minnesota potential.
The four points for each $u$ value correspond to calculations with
$K_{\rm max}=6,8,10$ and $12$, from the lowest to the largest value of $\alpha_E$,
respectively.
For $K_{\rm max}\geq 8$ we clearly see a linear dependence between
  $\alpha_E$ and $r_{\rm skin}$ for all the three $u$ values
where the coefficients depend on the separation energy as
\begin{equation}
\label{linear}
\alpha_E=a(S_{2n})+b(S_{2n})r_{\rm skin}\,.
\end{equation}
Because  $S_{2n}$ is converged and because of the linear dependence displayed
in Fig.~\ref{fig_alpha_skin} we deduce that the
calculation of $\alpha_E$ is
not fully converged because the radii, and especially $r_n$, are not fully
converged. The calculation of a radius of the ground state does not require an
expansion on the dipole excited states as in Eq.~(\ref{start}) and as such is
less computation demanding and can be performed for larger model spaces, where
radii are better converged.
Thus, the approach we take to estimate $\alpha_E$ from our calculations is to
fit the coefficients $a(S_{2n})$ and $b(S_{2n})$ from the
$\alpha_E$ results in the available model spaces and, assuming that this
physical linear dependence  will be unchanged in larger model spaces, we will
use the coefficients to obtain $\alpha_E$ from a bound state calculation of
$r_{\rm skin}$.
Starting with $u=1.05$, so that  $S_{2n}$ is close to
  experiment, we fit the parameters $a$ and $b$ to the
  results of our calculations using the available values of $K_{\rm max}\geq 6$.
  We test this procedure on the
  available model space by varying the largest $K_{\rm max}$.
  For a model space with largest $K_{\rm max}=10$, we  obtain $a=-0.7\pm0.2$~fm$^3$
and $b=1.83\pm0.3$~fm$^2$ fitting to
three points, $K_{\rm max}=6,8,10$. Using these values and the value $r_{\rm
  skin}=$0.776 fm, calculated in the
next largest model space $K_{\rm max}=12$, our linear ansatz of Eq.~(\ref{linear})
 yields $\alpha_E=0.7\pm0.3$ fm$^3$.
 The calculated value of $\alpha_E$ from the hyperspherical harmonics
 expansion up to $K_{\rm max}=12$ is
 0.754 fm$^3$, which is within our estimated error band.
Now we will repeat this procedure utilizing our best three values
$K_{\rm max}=8,10,12$ (we omitted the $K_{\rm max}=6$ point as it does not fall
in line with the other points). The resulting values are
 $a=-1.27\pm0.04$ fm$^3$ and $b=2.62\pm0.05$ fm$^2$.
We then calculate $r_{\rm skin}$ up to the largest grand-angular momentum value
  accessible with our computational facility,
$K_{\rm max}=16$.
Using the corresponding $r_{\rm skin}=$0.82 fm and propagating the
fit errors on $a$ and $b$ in the linear ansatz, we obtain $\alpha_E=0.88\pm0.06$
fm$^3$. Concerning the skin radius,
one can clearly see from
Fig.~\ref{rskin_fit}, that convergence is approached.
Extrapolating these points with an exponential {\it ansatz} of the form~ $r_{\rm skin}(K_{\rm max})=r_{\rm skin}({\infty})-c e^{-\kappa K_{\rm max}}$
we get $r_{\rm skin}(\infty)=0.87(5)$ fm. As an error estimate
we take the difference between $r_{\rm skin}(K_{\rm max}=16)$ and  $r_{\rm skin}(K_{\rm max}=12)$.
The theoretical value is somehow larger than the experimental data.
In fact, combining different measurements of the matter radius \cite{Tan92,Alk97,Kis05} with the most recent evaluation  of the proton radius \cite{Maxime},  one can infer $r_n$ and consequently the skin radius, which is found to be  0.52$\le r_{\rm skin}^{\rm exp}\le$0.62 fm.
The variation on $r_{\rm skin}^{\rm exp}$ is fairly large, due to the large
uncertainty in the matter radius determination from ion scattering.

Because the extrapolated skin radius is our best estimate of this observable, we use this value in
Eq.~(\ref{linear}) and propagate its error
estimating the polarizability, considering it independent from the fit errors on $a$ and $b$.
Finally, our estimate of the theoretical nuclear electric polarizability of $^6$He is
$\alpha_E=1.00(14)$ fm$^3$. This value is consistent  with
what we obtained without extrapolating the radius, showing that the error bars are based on conservative estimates.
 If we were to use Eq.~(\ref{linear}) with the
experimental values of the skin radius, one would obtain a nuclear electric polarizability of  $\alpha_E=0.08-0.34$ fm$^3$, which is even smaller than the estimate based solely on theory.

\section{Comparison with experiment}

\begin{figure}
\includegraphics[scale=0.35,clip=]{integrand.eps}
\caption{(Color online) Double inverse energy weighted cross section as a function of the energy $\omega$ and sum rule $m_{-2}(\bar{\omega})/2\pi^2$ as a function of $\bar{\omega}$ for $^6$He with the Minnesota potential and $u=1.05$.
}
\label{fig_integral}
\end{figure}

\vspace{0.5cm}
\begin{figure}
\includegraphics[scale=0.35,clip=]{6He_4He_comparison.eps}
\caption{(Color online) Double inverse energy weighted cross section as a function of the energy $\omega$ for $^4$He and $^6$He with the Minnesota potential and $u=1.05$.
}
\label{fig_comp}
\end{figure}

In Ref.~\cite{Moro} an evaluation of the experimental number for the
polarizability of $^6$He was presented, leading to  $\alpha^{\rm exp}_E=1.99(40)$ fm$^3$.
This was obtained from the inverse energy weighted integral of experimental and  theoretical $B(E1)$ response functions for $^6$He. The experimental distribution was measured by Aumann {\it et al.}
\cite{Aumann} from the Coulomb breakup of $^6$He on Lead and Carbon at 240 MeV/$u$ up to 8 MeV above threshold. In order to obtain the polarizability, data were extrapolated up to 12.3 MeV, where the threshold for the breakup into two tritons opens up. The theoretical curve was taken from a calculation of the dipole transition to the $1^-$ continuum  \cite{Thompson}
 in a three-body model with phenomenological $n-n$  and $n-\alpha$ interactions plus an effective three-body force, also calculated up to the two tritons threshold.
The estimate was done in two steps: ({\it i}) an average of the experimental data and theoretical curve was taken up to 12.3 MeV; ({\it ii})
to account for the higher energies, the polarizability of $^4$He was added. The latter one basically comes from integrating the photo-dissociation data  from Arkatov {\it et al.} \cite{Arkatov}.
Because we can access the full response functions using the LIT method at any energy below the pion production threshold, we can verify these two approximations. First, it is interesting to study
 the convergence of $\alpha_E$ calculated as a sum rule of the response, see Eq.~(\ref{moments}) and (\ref{alphaEPSR}), to investigate the validity of ({\it i}).
In  Fig.~\ref{fig_integral},  we present both the integrand function $\sigma_{\gamma}(\omega)/2\pi^2\omega^2$ versus $\omega$ and the convergence of the integral $m_{-2}(\bar{\omega})/2\pi^2$ versus $\bar{\omega}$.
At $\bar{\omega}=8$ MeV the sum rule is exhausted only up to 75$\%$. Thus only 75$\%$ of the $\alpha_E^{\rm exp}$ is based solely on experimental data.
At $\bar{\omega}=$12.3 MeV, where the two $^3$H channel opens up  the sum rule is exhausted up to
90$\%$. We observe that one needs to integrate up to $40$ MeV of energy to have the sum rule exhausted at the $98\%$ level.

To verify the approximation ({\it ii}) we can compare the
integrand function $\sigma_{\gamma}(\omega)/2\pi^2\omega^2$ for $^6$He and $^4$He at energies beyond 12.3 MeV. In Fig.~(\ref{fig_comp}) we observe that the two curves agree with each other for $\omega>35$ MeV, where the sum rule is almost exhausted. In the region beyond the $^4$He disintegration threshold and below about $35$ MeV one would over estimate the sum rule integrating the $^4$He curve, because  one gets 0.044 fm$^3$, to be compared to the 0.027 fm$^3$ obtained when correctly integrating the $^6$He curve. On the other hand, neglecting the part of the cross section for $\omega>12.3$ MeV and below the $^4$He disintegration threshold one underestimates the sum rule. The  contribution of this portion is  0.037 fm$^3$, about 5$\%$ of the sum rule. These two effects almost cancel out so that, the approximation ({\it ii}) does not lead to a big error.

We think that the main reason of the disagreement between the estimate from Ref.~\cite{Moro} and our
calculations comes from the difference in the low-energy part of the
response. In our previous work \cite{He6Sonia1, He6Sonia2} we have shown that
our calculations with semi-realistic potentials underestimates the data from
Aumann {\it et al.}~\cite{Aumann}. Thus, what we observe for the polarizability is consistent
with this fact. We would like to point out that (i) nuclear corrections might
affect the results in the ion scattering experiment of \cite{Aumann} and that (ii) as
discussed earlier, similar experiments for $^3$He
lead to large discrepancy with photo-dissociation results.
Nevertheless, measuring $\alpha_E$ from the dipole response
function it would be desirable to have data up to higher energies than 12.3
MeV.
Additional or alternative measurements of $\alpha_E$ would help to better constrain this observable.

\section{Conclusions}
We summarize our results as follows. We have carried out an estimate
 of the nuclear polarizability of $^6$He based on the hyperspherical harmonics expansion
with simple semi-realistic potentials.
Our calculations clearly reproduce the fact that the polarizability of the halo nucleus of $^6$He is much larger than that of the tightly bound $^4$He.
For $^4$He the semi-realistic Minnesota potentials lead to a value of the polarizability which is consistent with realistic calculations and is about 15$\%$ smaller than the average value in the experimental band. Nevertheless, a large disagreement is found for $^6$He.
In order to estimate $\alpha_E$
we have chosen a potential that reproduces the separation energy and then we investigated the correlation  of the polarizability with the skin radius.
Our final result is
$\alpha_E=1.00(14)$ fm$^3$,
which is about a factor of 2
smaller than the estimates from experimental data. This points towards a disagreement of microscopic theory and
experiments. To shed light on this, it would be nice to have more data or
alternative measurements of $\alpha_E$. Concerning the theoretical
calculations, it is desirable
to extend these results to realistic potentials including also three-body forces. We leave this subject to a future work.
\vspace{0.5cm}

\section{Acknowledgment}

The work of  R.~Goerke and S.~Bacca
was supported in part by the
Natural Sciences and Engineering Research Council (NSERC) and in part by the
National Research Council of Canada.
The work of N. Barnea was supported by the ISRAEL SCIENCE FOUNDATION
(grant no.~954/09).
Numerical calculations were  performed at TRIUMF.


\section{Introduction}
  \input{section1/intro}

\section{GNURadio code}
  \input{section2/tunnel}

\section{Results and future work}
  \input{section3/results}

\bibliographystyle{IEEEtran}


\subsection*{Results}

Here, we overcome these obstacles by mixing small amounts of nanodiamonds from solution directly into a known
photoresist \cite{Fischer2011} (See Methods). This special photoresist allows for the fabrication
of transparent optical elements made from an acrylate polymer containing nanodiamonds. The nanodiamonds contain
NV-centres which are photostable even after exposure by the DLW laser. After the DLW process they serve as
integrated single-photon sources. It is especially noteworthy that our technique is not limited to NV-centres in
diamond nanocrystals. It can be applied to any photostable single-photon emitter in nanocrystalline material,
e.g. other diamond defect centres, such as Silicon vacancy-centres \cite{Neu2011}, or other deep defect centres
 in large band gap material \cite{Weber2010}.

In the following, we demonstrate the fabrication, functionality, and interaction of key building blocks for quantum photonic
circuitry. These are interconnects (i.e., waveguides) and functional elements like resonators and emitters. We
measure single-photon emission from a diamond nanocrystal in a waveguide-coupled resonator which demonstrates
the functional combination of all these building blocks and paves the way to more complex devices.

As a first component for integrated photonic circuits, disc resonators with a disc diameter of
\unit{20}{\micro\metre} and a disc thickness of approximately \unit{1.2}{\micro\metre} on a stem with diameter
 \unit{10}{\micro\metre} were fabricated using the photoresist functionalised with nanocrystals (see Methods).
The resonator supports whispering gallery modes which are confined in the
disc close to the rim with high quality factors (Q-factors).
We measured the Q-factor using a tunable laser coupled into the resonator via a tapered
fibre\cite{Vahala2003}. We found Q-factors as large as $10^4$ at a wavelength of around \unit{770}{\nano\metre} (see Figure
\ref{fig:tapercoupling} and Methods). The observed free spectral range of \unit{6.5}{\nano\meter} of the
resonator corresponds well to the expected value of \unit{6.4}{\nano\meter} for a disc with diameter \unit{20}{\micro\meter}
derived from geometrical considerations assuming an index of refraction of 1.5.

\begin{figure*}
  \includegraphics{Figure_2}
  \caption{\textbf{Mode measurements}. \textbf{a}, A tunable laser is coupled to a disc resonator via a tapered fibre. Upon sweeping the
  laser frequency, the transmitted light is modulated by the modes of the resonator. \textbf{b}, Resonator mode with a quality factor Q of
  $1.2 \times 10^4$. \textbf{c}, Scan over many modes with a free spectral range of \unit{6.5}{\nano\meter}
  indicated by the red arrows.}
  \label{fig:tapercoupling}
\end{figure*}

Next, a homebuilt confocal microscope (see Methods) was used to raster scan the fabricated resonator discs  (see Figure 3a) and
to identify single NV-centres by measuring the second-order autocorrelation function $g^{(2)}(\tau)$ in a
Hanbury Brown and Twiss (HBT) interferometer. In parallel, a grating spectrometer was used to resolve the
emission spectra of the individual emitters. Figure \ref{fig:diainresonator}b shows a confocal scan of a
resonator where fluorescent defects can be identified as bright spots on the resonator. Encircled is a spot on
the resonator's outer rim where coupling to the disc's whispering gallery modes is expected. Figure
\ref{fig:diainresonator}d displays an autocorrelation measurement of fluorescence collected from that spot. A
clear antibunching dip is observed reaching a value of $g^{(2)}(0)=0.31\pm0.04$ (standard error) as deduced from the fit shown as red
curve. No background correction was applied to these data. This measurement shows that the bright spot indeed
corresponds to a single NV-centre in the disc resonator.

The fluorescence spectrum of an NV-centre at room temperature is typically broadened over \unit{200}{\nano\meter} from approx.
\unit{600}{\nano\meter} to \unit{800}{\nano\meter} by phonon sidebands. This corresponds well to the measured emission spectrum shown in Figure
\ref{fig:diainresonator}c. The peak at approximately \unit{630}{\nano\meter} stems from fluorescence of the
photoresist and bleaches after long excitation. For confocal scans and correlation measurements this background
light was suppressed by spectral filtering. In addition, there is a fine modulation of the spectrum due to the
cavity resonances. It is attributed to photons that were initially emitted into the resonator modes and
afterwards scattered out by the diamond nanocrystal. Hence, these photons are detected in addition to the flat
unmodulated spectrum emitted directly out of the disc.

\begin{figure*}
  \includegraphics{Figure_3}
  \caption{\textbf{DLW resonator containing single NV-centres.} A diamond nanocrystal containing a single NV-centre is coupled to whispering gallery modes
  of a DLW disc resonator. \textbf{a}, Measurement scheme. Detection and excitation take place at the same point in a confocal configuration.
  \textbf{b}, Scanning confocal image of the resonator disc. The circle indicates a bright spot identified as single NV-centre. Its
  fluorescence is analysed in \textbf{c} and \textbf{d}. Scale bar is \unit{5}{\micro\metre}. \textbf{c},
Spectrum of the collected fluorescence. The resonator modes are seen on the broad NV-centre phonon sidebands. The
peak at \unit{630}{\nano\meter} stems from the photoresist and can be bleached over time.
  \textbf{d}, Autocorrelation function of the fluorescence from the NV-centre. A clear antibunching behaviour can be seen. The red curve is a fit to the data.}
  \label{fig:diainresonator}
\end{figure*}

As a second photonic element, arc waveguides  from the same material with a width of \unit{1.8}{\micro\meter} and
a length of \unit{40}{\micro\meter} were built. The waveguides were fabricated close to the disc resonator
in order to enable evanescent coupling to the resonator's modes (see Figure \ref{fig:waveguide}a). Their ability
to guide light within a broad spectral range is demonstrated by illuminating one end of the waveguide with light from a halogen lamp as seen
in Figures \ref{fig:waveguide}b and c. The blue and red circular spots correspond to the input and output of
the arc waveguide, respectively.

Coupling via the waveguides to the resonator modes was confirmed by shining an excitation laser (wavelength of
\unit{532}{\nano\metre}) into one port of the waveguide and, after blocking excitation light by a dichroic
mirror, analysing the light coming out of the other port. Obviously, the \unit{532}{\nano\metre}-laser excites
background fluorescence in the resonator, which is then coupled back into the waveguide. The measured
fluorescence spectrum in Figure \ref{fig:waveguide}e thus shows the characteristic modulation by resonator
modes and proves a good coupling between guided modes in the arc waveguide and confined modes in the resonator
(see Supplementary Information).

\begin{figure*}
  \includegraphics{Figure_4}
  \caption{\textbf{Arc waveguide coupled to resonator.}
  \textbf{a}, Scanning electron micrograph of a 3D arc waveguide/resonator assembly. Length of the scale bar is \unit{5}{\micro\metre}.
  \textbf{b}, The waveguide resonator assembly as seen through the cover slip. The upper half of the image is illuminated from below by white light.
  Efficient light guiding through the arc waveguide can be seen. Scale bar is \unit{10}{\micro\metre}. \textbf{c}, Sketch of the
  region shown in \textbf{b} with the illuminated area marked yellow. \textbf{d}, Sketch how the spectrum in \textbf{e} was acquired.
  \textbf{e}, Spectrum of the fluorescence detected at one of the waveguide's ports while exciting through the other.
  Resonator modes can be clearly identified. \textbf{f}, Experimental configuration for the autocorrelation measurement in \textbf{g}.
  A second objective is used for excitation.
  \textbf{g}, Measured autocorrelation function of a NV-centre inside a waveguide. It is acquired at both output ports of the waveguide
  in a cross-correlation configuration.}
  \label{fig:waveguide}
\end{figure*}

By moving the excitation spot of a second objective over a waveguide (see Figure \ref{fig:waveguide}a), individual nanodiamonds inside
can be addressed and identified as spots with higher fluorescence. The intensity cross correlation
between the waveguide's two outputs shows a clear antibunching behaviour. This result highlights a way of
coupling single emitters to thin waveguides, which can possess very high collection efficiencies,
as demonstrated recently using  tapered optical fibres
\cite{Fujiwara2011,Schroder2012}. Here, we show a direct on-chip integration of such devices with both waveguide
ends accessible. Detection at the waveguide's ends can yield count rates of
 \unit{29,000}{\per \second} after passing a 50/50 beam splitter and spatial and spectral
filtering. The autocorrelation function at zero time delay of the fluorescence collected in this way is
$g^{(2)}(0)=0.37\pm0.13$ (standard error) showing the single-photon character of the light (see Figure \ref{fig:waveguide}g). When
accounting for the 50/50 beam splitter, at high excitation powers the number of background corrected counts
from the defect centre is measured to be \unit{130,000}{\per \second}.

\begin{figure*}
  \includegraphics{Figure_5}
  \caption{\textbf{Waveguide coupling of a single NV-centre inside a resonator.} \textbf{a}, Sketch of the experimental configuration.
  The excitation spot is scanned over the resonator disc. Photons are detected at both waveguide outputs simultaneously.
  \textbf{b}, Photon counts collected at
  one end of the waveguide while scanning the excitation spot with a second objective. The circle
  highlights the position of a single NV-centre. Shape distortions are due to non closed-loop piezo-scanning. Scale bar is \unit{5}{\micro\metre}.
  \textbf{c}, Cross correlation (background corrected) between the waveguide output ports when the NV-centre marked by a
  circle in \textbf{b} is excited. The bin size corresponds to the repetition rate of \unit{40}{\mega\hertz}}
  \label{fig:res16dia}
\end{figure*}

Finally, we verified active quantum functionality of a coupled 3D system, i.e., single photon emission to
resonator modes, as well as subsequent collection and routing via a coupled waveguide. In order to do so, raster
scans of the excitation spot over the resonator disc were performed while the photons from one of the
waveguide's ends are collected (see Figure \ref{fig:res16dia}a). When the excitation laser hits the outer rim of
the resonator, the created fluorescence light is coupled to the resonator modes and subsequently out to the
waveguide, as can bee seen in Figure \ref{fig:res16dia}b. Therefore, the outer rim of the resonator appears
bright in the raster image. Bright spots indicate possible positions of NV-centres. To identify NV-centres via
their non-classical light emission, an HBT setup, with the laser-written arc waveguide serving as beam splitter, is
used.
Figure \ref{fig:res16dia}b shows this correlation measurement for the NV-centre encircled in Figure
\ref{fig:res16dia}a. The data shown were acquired with a pulsed laser at \unit{40}{\mega\hertz} repetition rate and background
corrected (see Supplementary Information). The value of $g^{(2)}(0)=-0.18\pm0.21$ (standard deviation) is far below the threshold
$g^{(2)}(0)=0.5$ for having
equal contributions of two centers and proves that the main contribution
to the signal stems from a single NV-centre. This clearly indicates that the NV-centre emits single photons
into resonator modes, which are then coupled out and redirected by the arc waveguide. With a single source in a resonator
coupled to a waveguide, this device represents a key integrated 3D quantum photonic circuit.

\subsection*{Discussion}

In conclusion, we have shown an approach which allows for the direct incorporation of single quantum emitters
into true 3D photonic structures of nearly arbitrary shape. The one step fabrication scheme is easy to implement
and requires no elaborate clean-room environment. We have pointed out that our method is not limited to
NV-centres in nanodiamond, but can be extended to any other stable quantum emitter embedded in nanocrystals. The
flexibility of DLW enables immediate scaling-up to more complex structures (see Figure \ref{fig:dlw}b) where
single photons are collected, sent through beam splitters, interferometers, or other optical elements
inter-connected by waveguides in a 3D architecture.
Also tree-dimensional multi-arm quantum interferometers~\cite{} can be implemented.
Tuning of individual components, e.g. the resonator, is
possible as well (see Supplementary Information). With (partial) metallisation of the structures
embedding single photon emitters in 3D plasmonic structures or in metamaterials \cite{Rill2008} can be envisioned.
A next challenging, but rewarding step is to perform DLW of perfectly aligned 3D structures around emitters of known position.
This would enable on-demand fabrication of arbitrary 3D quantum photonic architectures.

\subsection*{Methods}

\textbf{Direct laser writing.} The DLW setup consists of a Ti:sapphire oscillator (SpectraPhysics MaiTai HP)
delivering \unit{100}{\femto\second} pulses centred around \unit{810}{\nano\metre}. The beam is focussed through
an oil-immersion lens (Leica HCX PL APO 100x/1.40-0.70 OIL CS). The sample is moved with a 3D piezo stage (PI
P-527.3CL). The writing was performed with \unit{6}{\milli\watt} average writing power and at
\unit{50}{\micro\metre\per\second} scan velocity. The samples were developed in isopropyl alcohol for 30
minutes, afterwards in acetone for 30 minutes, and finally dried super-critically (Leica CPD030).

The new photoresist is based on the monomer pentaerythritol tetraacrylate (PETTA) which contains 350~ppm
monomethyl ether hydroquinone as inhibitor. Next, 0.25~\% wt of the photoinitiator
7-diethylamino-3-thenoylcoumarin were added. Last, 2~\% wt of an ethanol-based nanodiamond suspension were added
and stirred overnight. Nanodiamonds were of type 1b and had a median diameter of \unit{25}{\nano\metre} (Microdiamant AG).
Spatial dimensions of the fabricated structures were extracted from scanning electron micrographs.

\textbf{Confocal microscope.} The homebuilt confocal microscope is based on the body of a Zeiss Axiovert 200.
Either a frequency doubled Nd:YAG continuous wave laser or a pulsed laser (LDH-P-FA-530, PicoQuant) were used to
create the excitation light around \unit{532}{\nano\metre}. A large-numerical-aperture air microscope objective
mounted on a 1D piezo stage (MIPOS 100, piezosystem jena) focussed the excitation light and collected the
fluorescence. The sample was scanned by a 2D piezo stage (PXY 80 D12, piezosystem jena). Collected fluorescence
light was separated from the excitation laser by a dichroic mirror. Then, it was either spatially filtered and
sent
to a grating spectrometer
or it was spectrally filtered by RG630 glass, split by a 50/50 beam splitter, and directed towards two pinholes.
These pinholes could be independently adjusted to spatially filter out different spots of the image plane, so
that cross correlations could be measured. After the pinholes the light was directed to avalanche photodiodes
(SPCM-AQR-14/SPCM-AQRH-14, Perkin Elmer).

To perform correlation measurements at the waveguide ends, an oil-immersion microscope objective (UPlanSApo
60XO, Olympus) was used for collecting the light, while different air objectives were used to focus the
excitation laser onto the sample from the other side. The air objectives were mounted on a 3D piezo system
(PXY100 ID and MIPOS 100, piezosystem jena) to enable scanning of the excitation spot while the sample was kept
at a fixed position. Also, the 50/50 beam splitter within the HBT setup could be replaced by a D-shaped mirror
in an image plane in cases where cross-correlation measurements with high count rates were needed.

\textbf{Q-factor measurements.} The mode structure of the disc resonators was analysed by means of coupling in a
tunable external-cavity diode laser (Velocity Series, New Focus) at around \unit{770}{\nano\metre} via the
evanescent fields of a tapered optical fibre.
By tuning the frequency of the laser over the distinct whispering gallery modes, light was coupled into the
resonator and different modes could be observed as Lorentzian shaped dips in the transmitted power. The
polarization of the incoming light was chosen to maximize the coupling depth. For normalizing the data sets, we
first performed a reference scan with an uncoupled fibre taper and compared the results to the situation when
the same taper was coupled to the resonator. This provided direct access to the different coupling depths and
cancelled out overlaying power modulations caused by slightly non-adiabatic taper transitions.

The Q-factors were calculated from the dips by a Lorentzian fit function with an additional linear term to
better match the local environment of the resonances. The fibre taper used had a waist diameter of
\unit{1.5}{\micro\metre}. Measurements were performed in full contact.

\subsection*{Acknowledgements}
The Berlin team acknowledges support by DFG (FOR1493 and SFB951).
The Karlsruhe team acknowledges support by the
DFG-Center for Functional Nanostructures (CFN) via
subprojects A1.4 and A1.5 and by the Karlsruhe School of Optics \& Photonics (KSOP).
J.W. acknowledges funding by the state of Berlin (Elsa-Neumann).

\subsection*{Author contributions}
A.W.S. wrote the manuscript and conceived the experiment with J.F., J.K. and J.F. produced the samples, A.W.S., R.H. and J.W.
performed the optical measurements and data analysis. M.W. and O.B. supervised the study.
All authors discussed the experiment, the results and the manuscript.

\subsection*{Additional Information}
The authors declare no competing financial interests.