text:
\section{ Introduction}

The partition functions for conformal field theories  on surfaces with boundaries are
of interest in diverse areas of physics.
   In soft condensed matter physics it is of interest to calculate partition
functions on the Euclidean  plane  with compact boundaries for the purpose of understanding fluctuation induced forces \cite{kardar}.
 In string theories, amplitudes correspond  to such partition functions  and one must sum over all moduli
 and Euler characteristics. However, to calculate  a particular gauge slice
is not always simple. Powerful methodologies in boundary conformal field theory (BCFT)
 \cite{Cardy1} have been utilized to determine the spectrum of primary operators in Gaussian as well as non-Gaussian   theories.
 These methods utilize the state operator correspondence to implement the boundary conditions
 via the insertion of  boundary operators.

In this paper  we will take a different approach based on world-line effective field theories (EFT)\cite{GnR} as applied
to statistical systems \cite{YRD1,eisenriegler}.
In particular,  we will be interested in the case of partition functions on the Euclidean plane with
compact boundaries. The EFT approach is similar to  BCFT  in that
the effects of the boundaries are accounted for via operator insertions.  However, we will make no use of
the state operators correspondence, and as such, the technique will not be limited to conformal field theories.

We will begin by considering
the free bosonic action
\begin{equation}
S_0^B =\frac{1}{2} \, \! \int d^2r\,\partial\phi
 \partial \phi
  \ , \label{H}
\end{equation}
on a plane with holes which satisfy
Neumann boundary conditions. This problem is identical to finding the partition function for
a film embedded with rigid disks which are allowed to bob but  not tilt\footnote{This is not true for all observables but will hold for the partition source-free partition function.},
where the notion of tilting is given meaning by a particular choice of
gauge. In particular the scalar field partition function will correspond to
the Monge gauge choice \cite{CMTtext} where the height of the field $\phi$ is measured with respect to a Euclidean flat base plane.  In the film model there are non-linear corrections
which are there as a consequence of reparameterization invariance, but it can
be shown that these corrections are suppressed \cite{YRD1} by powers of $k/ \beta$
where $k$ is the modulus characterizing the surface energy.
The fluctuation induced forces  in these systems are  relevant to soft condensed
matter problems \cite{kardar}.

In \cite{YRD1} it was shown that it is possible to utilize the world-line effective field theory
developed for studying black hole inspirals \cite{GnR} to study the Casimir forces on two dimensional surfaces (the dimensionality is actually irrelevant).  The idea is to work in the point particle
approximation, and then account for the finite size (read boundary conditions) by
adding an infinite tower of higher dimensional operators to the action.
For defects with only linear responses to external perturbations   these terms are
all quadratic in the field, and their dimensionality is fixed by the number of
derivatives. Power counting shows that non-linear interactions on the world-line
generate corrections which are suppressed at high temperatures.
 In  the limit where the size of the
objects ($R$) is small compared to the distance between them $r$, truncating the
infinite sum induces errors of order $(R/r)^\tau$, where $\tau$ is the the number of
derivatives in the first term dropped in the series.
The coefficients of these operators (Wilson coefficients)  correspond
to generalized polarizabilities and are fixed by performing a matching calculation
between the full theory and the effective theory.  The full set of Wilson coefficients can be calculated
 by first matching the response of  the individual
 holes to multipole fields in the full and effective theory. We will then show that these matching
 coefficients can equally well be extracted by determining the boundary
 operators in the Hamiltonian version of the theory.
Comparing the two procedures
 illuminates a map between closed string excitations and effective polarizabilities.
The model is then solved by summing the series for the higher dimensional
operators into a non-local operator which can be used to
write down a closed form expression for the partition function.

In the effective theory the action for  for n-bodies may be written as
\begin{equation}
S=S_0+\frac{1}{2}\sum_{i} C_i (\partial_{a_1}....\partial_{a_n} \phi)(\partial_{a_1}....\partial_{a_n} \phi),
\end{equation}
where the second term is localized at the origin.
Operators involving traces over indices are redundant in this theory as they
can be traded for other operators via a canonical transformation.
For non-spherically symmetric objects the polarizibilities will transform
non-trivially under the two dimensional rotation group, but are singlets
for the spherical  case of interest here. These coefficients are fixed by a simple
matching procedure whereby we place the objects in a background
field and calculate its response.  For a given polarizability of a given
order
we may extract the corresponding $C_i$ by placing the object in a multipole
field of the same order.

To utilize holomorphy we match in the complex plane a general n-pole  background
field takes the form
\begin{equation}
\label{bg}
\phi_{bg}(z,\bar z)= \frac{\lambda}{R^n} (z^n+\bar z^n).
\end{equation}
We will impose Neumann boundary conditions which
preserve the conformal symmetry.
\begin{equation}
\partial_r \phi(z)\mid_{R}=0
\end{equation}
  This choice means that
the first non-zero $C_i$ will be $C_2$, a dipole polarizability.
In the full theory we solve the rudimentary boundary value problem.
The induced field is given by
\begin{equation}
\delta\phi= \lambda  R^{n}(\frac{1}{z^n} +\frac{1}{\bar z^n}).
\end{equation}

Expanding around the background field (\ref{bg})
utilizing the Greens function in complex coordinates
\begin{equation}
G(z,\bar z) = \frac{1}{4 \pi } (\log (z) +\log (\bar z))
\end{equation}
and calculating the one point
function leads to the EFT result
\begin{equation}
\delta \phi =-\frac{C_n \lambda}{4 \pi  R^{n} } 2^{2n-1} n! (n-1)!(\frac{1}{z^{n}}+\frac{1}{\bar z ^{n}}).
\end{equation}
We may then extract the matching coefficient to be
\begin{equation}
\label{Cs}
C_n= -\frac{\pi }{n! (n-1)!~2^{n-2}}R^{2n}.
\end{equation}

This form of matching coefficient can also be obtained  from perspective of closed string
theory. If we conformally map an isolated hole  into an infinite cylinder, and consider
the direction along the length as being time-like, then  the partition function corresponds  to a closed string propagator
emanating from a  boundary state \cite{callan}.
The boundary states are annihilated by the constraint
\begin{equation}
\partial_1 \phi \mid \psi \rangle =0
\end{equation}
where the world-sheet is parameterized by $(\sigma_1,\sigma_2)$ and $\sigma_1$ is the
time like direction.

Our gauge choice (Monge) corresponds to choosing mode oscillators in the string expansion to be  along one space like direction in the embedding space.
 In terms of the mode operators the constraint \cite{callan}reads
\begin{equation}
(\alpha_n+\tilde  \alpha_{-n} )\mid \psi \rangle =0
\end{equation}
where the (standard) normalization is fixed via
\begin{eqnarray}
\partial \phi(z) &=& -i \left( \frac{\alpha^\prime}{2}\right)^{1/2}\sum_{m=-\infty}^{\infty} \frac{\alpha_m}{z^{m+1}},\end{eqnarray}
\begin{equation}
[\alpha_m,\alpha_n]=[\tilde \alpha_m,\tilde \alpha_n]=m \delta_{m,-n},
\end{equation}
and
\begin{equation}
S=\int \frac{1}{4\pi \alpha^\prime} (\partial \phi \bar \partial \phi) d^2 \sigma.
\end{equation}

The solution to this constraint is
\begin{equation}
\mid \psi \rangle  \propto  e^{-\sum_{p=1}^\infty (\frac{1}{p} \alpha_{-p} \tilde \alpha_{-p})}\mid 0 \rangle.
\end{equation}
Using the fact that
\begin{eqnarray}
\alpha_{-m} &=& \left( \frac{2}{\alpha^\prime}\right)^{1/2} \frac{i}{(m-1)!}\partial^m \phi(0)
\end{eqnarray}
the state is written as
\begin{equation}
\mid \psi \rangle  \propto  e^{-\sum_{p=1}^\infty \frac{2}{\alpha^\prime p! (p-1)!}(\partial^p\phi \bar \partial^p \phi)}\mid 0 \rangle.
\end{equation}
rescaling the fields so that they are normalized according to (\ref{H}) leads to
\begin{equation}
\mid \psi \rangle  \propto Exp\left[{-\sum_{p=1}^\infty \frac{4\pi  R^{2p}}{ p! (p-1)!}(\partial^p\phi \bar \partial^p \phi)}\right]\mid 0 \rangle.
\end{equation}
Taking this as an operator insertion
in the path integral
reproduces the Wilson coefficients $C_n$. It is interesting to note that there is a
correspondence between the string excitation and the generalized susceptibilities.
For instance the dipole polarizability $(\partial \phi)^2$ on the world line maps to the dilaton excitation of the string. Higher order multipoles map to massive string excitations.

Given the action the two body problem can be solved.
We utilize the holomorphic nature of the fields to re-write the action  as
\begin{eqnarray}
S_{int}
&=& \sum_n C_n (2)^{n}(\partial_{z}^n \phi ) (\partial_{\bar z}^n \phi )\nn \\
\end{eqnarray}
Then using the result  (\ref{Cs}) we have
\begin{eqnarray}
&&S=\int d^2z  [\frac{1}{2}\partial_z \phi \partial_{\bar z} \phi+
\sum_a  (4 \pi )  R_a \sqrt { \partial_z \partial_{\bar z}}I_1(2R_a \sqrt { \partial_z \partial_{\bar z}})\phi(z) \phi(\bar z)\delta^2(z-z_a,\bar z - \bar z_a)].\nn \\
\end{eqnarray}
The exact result for the partition function follows from treating these operators as mass insertions and
recognizing that the contractions between two identical world lines lead to trivial
renormalizations that are pure counter-terms. Thus the only relevant vacuum diagrams are those
in which no two mass insertions are adjacent.

Let us consider the case of two particles.   To calculate the free energy we draw all possible vacuum diagrams
which have no disconnected sub-components.  We must always add two insertions at a time since
for  an odd number of insertions there would have to be a self-energy graph inserted which we renormalize to zero.
Diagrams which are related by a
reflection symmetry are not distinct.
\begin{figure}[thb]
\centering
\includegraphics[scale=.89]{fig1}
\caption{All connected vacuum diagrams can be written in term of this
topology. The different colored dots represent insertions at each of the two particles respectively.}
\label{fig1}
\end{figure}
At second order in the interaction the diagram reduces down to a simple bubble
and the resulting free energy is given by
\begin{eqnarray}
\label{exact}
-\beta F  &=&
\sum_{n=1}^{\infty}\prod_{a=1,3,5..}^{2n-1} (R_{A}R_{B})^n\frac{1}{n} \sqrt{ \partial_{z_a}\partial_{\bar z _a}}\sqrt{ \partial_{z_{a+1}}\partial_{\bar z _{a+1}}}I_1(2R_A\sqrt{ \partial_{z_a}\partial_{\bar z_a}})I_1(2R_B\sqrt{ \partial_{z_{a+1}}\partial_{\bar z_{a+1}}}) \nn \\
&\times& \log( z_a- z_{a+1})\log(\bar z_{a}-\bar z_{a+1})\nn \\
\end{eqnarray}
Each $n$ represents two insertions, one of each type, and $z_{2n+1}=z_1$.
Upon expanding this expression and differentiating one sets $z_a=z_1$ and $z_{a+1}=z_2$.
Let us now expand this result in $R/r$ to compare with the literature on
the lower order contributions.
The  leading order  contribution is the dipole-dipole interaction.
Expanding out (\ref{exact}) gives
\begin{eqnarray}
-\beta  F^{(2)}_{d-d} =
&=& R_A^2R_B^2  \partial_{z_1}\partial_{\bar z _1} \partial_{z_{2}}\partial_{\bar z _{2}} \log((z_1-z_{2})\log(\bar z_1-\bar z_{2}))\nn \\
&=&\frac{1}{r^4}R_A^2R_B^2
\end{eqnarray}

We also generate  non-linear terms such as the case  with two dipoles of each type in which case $i=2$.
\begin{eqnarray}
-\beta F^{(2)}_{d^2-d^2} = && R_{A}^2R_{B}^2\frac{1}{2} \sqrt{ \partial_{z_1}\partial_{\bar z _1}}\sqrt{ \partial_{z_{2}}\partial_{\bar z _{2}}}I_1(2R_A\sqrt{ \partial_{z_1}\partial_{\bar z_1}})I_1(2R_B\sqrt{ \partial_{z_{2}}\partial_{\bar z_{2}}}) \log(\bar z_1-\bar z_{2})\log(z_1-z_{2})\nn \\
&\times& \sqrt{ \partial_{z_3}\partial_{\bar z _3}}\sqrt{ \partial_{z_{4}}\partial_{\bar z _{4}}}I_1(2R_A\sqrt{ \partial_{z_3}\partial_{\bar z_3}})I_1(2R_B\sqrt{ \partial_{z_{4}}\partial_{\bar z_{4}}}) \log(\bar z_3-\bar z_{4})\log(z_3-z_{4})\nn \\
&=&
R_{A}^4R_{B}^4\frac{1}{2r^8} \nn \\
\end{eqnarray}
All of these perturbative results agree with those in the  literature \cite{GGK,LO}.
New results can be read off by going to higher order.

The case of N holes follows in a similar fashion. If we are interested in the N-body potential   we have a minimum of N insertions and we may add one
vertex insertion at a time.
\begin{eqnarray}
\label{sum}
-\beta F_N=&&\sum_{n=N}^{\infty}\sum_{\cal P}\prod_{a=1}^{n} R_1^{i_1} R_2^{i_2}.... R_N^{i_N} \sqrt{ \partial_{z_a}\partial_{\bar z _{a+1}}}I_1(2R_a\sqrt{ \partial_{z_a}\partial_{\bar z_{a+1}}})
 \times G( z_a- z_{a+1},\bar z_a- \bar z_{a+1})\nn \\
\end{eqnarray}
The second sum is over all possible choice non-equivalent permutation of  vertices.
Where two permutations are in the same equivalence class  if they can be rotated into each other.
 Defining $G(0)=0$ implies that all choices with two adjacent
identical vertices vanish. For each choice of vertices,  $(i_1,i_2,i_3...)$ correspond to the number of vertices of (hole) type ($1,2,3...$) in the
particular term in the sum. The expression (\ref{sum}) is easily evaluated to arbitrary order in MATHEMATICA, though
this particular form is most probably not the most efficient since it entails redundancies as compared to the two body
results presented above.

\section{CFT Solution for Two Holes}
The determinant of the Laplacian on Riemann surfaces  in many instances can recovered
by conformally mapping surfaces for which the partition is easily calculated to surfaces
which present more of a challenge, such as the case of interest in this paper.
The techniques used here are well known though it seems that the  explicit result for the case of interest, two spheres in a plane, is absent in the voluminous literature
on the subject of functional determinants in two dimensions.

The partition function on the plane with two holes can be calculated by a series of conformal
maps starting from the cylinder. One could of course map directly from the cylinder
to the plane with holes, but it is perhaps simpler to use the annulus as an intermediate step.
In each step we work in a conformal gauge
\begin{equation}
g_{ab} \sim e^{2 \phi} \hat g_{ab},
\end{equation}
in which case the following quantity is an invariant  \cite{alvarez,weisberger}
\begin{equation}
I_N
= \frac{Det^\prime_N(\nabla^2)}{A}( e^{\frac{L}{6\pi} -\frac{1}{4\pi} \int_{\partial M} \kappa ds}).\end{equation}
 Where we have chosen Neumann boundary conditions for the purpose of illustration.
 $L$ is defined as
\begin{equation}
L(\phi,\hat g)= \int_M d\hat \mu \hat K \phi +\int_{\partial M} d\hat s \hat \kappa \phi+\frac{1}{2} \int_M d\hat \mu \hat g^{ab}(\partial_ a \phi \partial_b \phi),
\end{equation}
and  hatted quantities are defined by the hatted metric, $K$ is the Gaussian curvature on $M$, $\kappa$ is the geodesic curvature on the boundary, and $A$ is the area of $M$. The prime above the determinant denotes the fact that
the zero mode is excluded.

Beginning with the partition function on a cylinder with length $l$ and unit radius. A rudimentary calculation yields the  result\begin{equation}
Z= \frac{1}{\sqrt{2l}} e^{\frac{\pi l}{6}} \prod_{n=1}^\infty \frac{1}{1- e^{-4\pi l n}}.
\end{equation}
We then map this result to the annulus using the conformal transformation
\begin{equation}
z= r_2 e^{2 \pi i w}
\end{equation}
with
\begin{equation}
2\pi l = \ln(r_2/r_1).
\end{equation}
The outer (inner) radii of the annulus are $r_2 (r_1)$. The conformal factor is given by
\begin{equation}
\phi= \frac{1}{2}\log( 2 \pi z \bar z).
\end{equation}
Given that both spaces are flat we find
\begin{equation}
L_{ann}= \pi \log (r_2/r_1).
\end{equation}
The annulus partition function is then given by \cite{weisberger}
\begin{equation}
Z_{ann} = \frac{f(r_1^2/r_2^2)}{\sqrt{2(r_2^2-r_1^2)}} (r_1/r_2)^{-1/6}
\end{equation}
\begin{equation}
f(x)= \prod_{n=1}^\infty\frac{1}{1-x^{n}}
\end{equation}
Then we can arrive at case of interest  (z plane with two holes)  via a conformal map from the annulus (w plane).

\begin{figure}[thb]
\centering
\includegraphics[scale=.93]{map}
\caption{Mapping between annulus and twice punctured plane}
\label{map}
\end{figure}
Tha mapping   shown in  figure (\ref{map})  is given by
\begin{equation}
\label{Cmap}
w=\frac{z-a}{a z-1}
\end{equation}

 where
 \begin{equation}
 a= \frac{1+x_1 x_2 +\sqrt{(x_1^2-1)(x_2^2-1)}}{x_1+x_2}
 \end{equation}
 \begin{equation}
 R= \frac{x_1 x_2 -1 -\sqrt{(x_1^2-1)(x_2^2-1)}}{x_1-x_2},
 \end{equation}
under the conditions
$x_2<a<x_1$.

 The metric on the  annulus is then mapped to the flat z plane
 \begin{equation}
  dz d\bar z=d w d\bar w \frac{(a^2-1)^2}{(aw-1)^2 (a \bar w-1)^2}
 \end{equation}
 so that
 \begin{equation}
 \phi=\frac{1}{2} \ln \left(  \frac{(a^2-1)^2}{(aw-1)^2 (a \bar w-1)^2} \right)
 \end{equation}
 and $\hat g=\delta_{ab}$ is the metric on the flat plane so that $\hat K=0$.
 The annulus is flat ($R \sim \nabla^2 \phi$) so $\int_{\partial M} \kappa =0$ since its Euler characteristic  $\chi=0$.
 However, the integral
 \begin{equation}
 \int_{\partial M} d\hat s \hat \kappa \phi=\int_0^{2\pi} d \theta (\phi(1,\theta)-\phi(R,\theta))= 2 \pi Log[R]
 \end{equation}
 does not vanish since $\phi$ is not constant along the boundary.

 The bulk contribution to the anomaly must vanish since the conformal map is an element of the $SL(2,C)$ subgroup
 with vanishing Schwarztian. Hence we are left the the following partition function for the plane with two holes.
 \begin{equation}
 \label{CFT}
 -\beta F= Log[ f[ R^2]]
 \end{equation}
 where
 \begin{equation}
f[x]\equiv  \prod_{n=1}^{\infty}\frac{1}{ (1- x^n)}
 \end{equation}
and
 \begin{equation}
 R= 1/2(r^2-r \sqrt{r^2-4} -2).
 \end{equation}
 Expanding in the ratio of the size of the holes relative to the distance between them gives
 \begin{equation}
 \beta F \approx \frac{1}{r^4} + \frac{4}{r^6}+\frac{31}{2r^8 }+\frac{60}{r^{10}}+...
 \end{equation}
 Which agrees with the expansion of (\ref{exact}).
 The expression (\ref{CFT}) certainly has many advantages over the result (\ref{exact}). In particular, though the
 expansion in (\ref{exact}) has a finite radius of convergence we can not expand it around to
 the osculating limit. Whereas  the Dedekind function which results from (\ref{CFT},41) has multiple
 representations which can be easily expanded around this limit. The result (\ref{CFT}) may be used
 be compared to
 the so-called ``proximity force approximation''\cite{NPA} which has recently been systematized in
 terms of a derivative expansion \cite{systematic}.
 \section{Non-Conformal Theories}
It is interesting to see that  the effective theory methodology can still
 yield exact results even for non-conformal kernels.  A non-conformally invariant  kernel of particular interest is the bi-Laplacian which
 corresponds to a model of a membrane where bending dominates tension.
  In this case there are more operators at each order in $R$
 since the equations of motion can no longer be used to eliminate operators involving a single Laplacian.
In the case of the
bi-Laplacian, beyond second order in derivatives
 we have two operators at each order in the derivative expansion.
At  the quadrapole  level we have two possible
operators.
 A convenient basis to choose is
\begin{equation}
O^{2}_a= \partial_z \partial_z \phi \partial_{\bar z} \partial_{\bar z} \phi,~~~O^{2}_b=  \partial_z \partial_{\bar z} \phi \partial_{\bar z} \partial_{ z} \phi.
\end{equation}

At higher orders we add one power of $\partial_{z,\bar z}$ to each side. So that at $n'th$  order in derivatives we have
\begin{equation}
O^{n+1}_a= \partial^{n+1}_z \phi  \partial_{\bar z}^{n+1}\phi,~~~O^{n+1}_b=  \partial_z^n \partial_{\bar z}\phi \partial_{\bar z} ^n\partial_{ z} \phi.
\end{equation}
In the case of the bi-Laplacian one must impose two boundary conditions to find unique
solutions.  Which operator has the first non-vanishing susceptibility will depend upon the particular choice.
The crucial point is that the number of operators at each level in derivatives does not grow which allows
again for  the matching procedure to carried out to all order in the derivative expansion \cite{YD}.

This fact may at first seem surprising given that the theory is not conformal. However,  the theory is invariant under the
$SL(2,C)$  sub-group of the Virasoro group  which is revealed once one makes a field redefinition \cite{Bludman,Ferrari}
\begin{equation}
\phi \rightarrow \mid\! \frac{dz}{dz^\prime} \!\mid \phi^\prime (z^\prime).
\end{equation}
In addition, higher order kernels can also be solved by the same method. However, as one
goes to higher orders more (but still finite) operators appear at each order.

\section{Conclusions}
Partition functions for on  higher genus surfaces  can have a rich structure.
The boundaries introduce new sets of operators which make even free theories
highly non-trivial.
 Models which are integrable on simple manifolds, such as
cylinders,  can be solved on more complicated surfaces through conformal mappings and
the associated anomalous transformations of the stress energy in the bulk and boundary \footnote{The partition function for sections of  spheres have been calculated in this fashion in \cite{dowker}}.
However, mappings to the Euclidean plane with multiple boundaries of various shapes
can become analytically intractable.

In this letter we have shown that two dimensional  partition functions
can also be expressed in terms of an infinite set of Wilson coefficients in a world line EFT. For Gaussian theories, such as the free bosonic one discussed here, these coefficients are exactly calculable. These coefficients which physically correspond
to generalized susceptibilities, are equal to the weights of string excitations in the coherent state
boundary operator in the closed string channel.
Given these coefficients, the partition function follows from the vacuum energy graphs. The result is
written as an infinite sum which is clearly inferior to  the infinite product
form which arise in the context of CFT methodology.
For the case of two spherical holes we have shown the equality of the CFT and EFT results expanding to arbitrary order in the separation
between the holes.

  However, as the number of
defects increases the only impediment to solving the theory using the world-line methodology is that the  series involved
grows more cumbersome. Here we only presented the sums in their most raw form,
and it is surely possible to simplify them further \cite{YRD}. Furthermore, it is also possible
to determine the Wilson coefficients for deformed holes using the fact that
can take advantage of rudimentary conformal transformation techniques to calculate
the infinite set of Wilson coefficients for ellipses  \cite{Haussman}.
We have also shown that  it is also possible to solve for the Wilson  coefficients for the  non-conformal bi-Laplacian kernel.

Finally it is interesting to ask how these EFT techniques can be applied to interacting theories.
For any conformal theory, the conformal mapping technique is still applicable, and there is
no reason that one could not find the fluctuation induced force say for an Ising model \cite{Eisenriegler}.
However, again if one is interested in more complicated geometries (such as multiple holes) the mappings will
become analytically intractable. However, given that distortions of single holes to other shapes, such as
ellipses, is fairly simple, the matching procedure in the effective theory technique is still relatively straightforward.
Thus, given that for a simple geometry one can solve for the Wilson coefficient, at least for an integrable
model, one can also determine the coefficients for more complicated shapes. The free energy for
multiple holes can then be calculated via the vacuum energy diagrams as shown above, at least for
 theories which are integrable.

\begin{acknowledgments}
This work  is supported by US DOE contract, 22645.1.1110173.
I thank the Caltech high energy group for its hospitality
and  the Gordon and Betty Moore Foundation for support.
I benefited from conversations with Anton Kapustin and Rich Holman, Markus Deserno and Cem Yolcu. I thank W. Goldberger for
comments on the manuscript and Erich Eisenriegler for pointing me to \cite{Eisenriegler}.
\end{acknowledgments}


text:
\section{Introduction}
Recently, observations of compact stars have provided new data of high
accuracy which put strong constraints on the high-density behaviour of the
equation of state of strongly interacting matter otherwise not accessible in
terrestrial laboratories.
In particular, the high masses of $M=1.96 +0.09/-0.12~M_\odot$ and
$M=2.73 \pm 0.25~M_\odot$ obtained in very recent measurements on the
millisecond pulsars  PSR B1516+02B and PSR J1748-2021B, repectively
\cite{Freire:2007jd}, together with the large radius of $R > 12$ km for the
isolated neutron star RX J1856.5-3754 (shorthand: RX J1856)
\cite{Trumper:2003we} point to a stiff equation of state at high densities.
Measurements of high masses are also reported for compact stars in low-mass
X-ray binaries (LMXBs) as, e.g.,
$M=2.0\pm 0.1~M_\odot$ for the compact object in 4U 1636-536
\cite{Barret:2005wd}.
For another LMXB, EXO 0748-676,
constraints for the mass $M\ge 2.10\pm 0.28~M_\odot$
{\it and} the radius $R \ge 13.8 \pm 0.18$ km
have been derived \cite{Ozel:2006km}.
The status of these data is, however, unclear since the observation of a
gravitational redshift $z=0.35$ in the X-ray burst spectra
\cite{Cottam:2002cu}
could not be confirmed thereafter despite numerous attempts
\cite{Cottam:2007cd}.

Measurements of rotation periods below $\sim 1$ ms as discussed for
XTE J1739-285 \cite{Kaaret:2006gr}, on the other hand, would disfavor too
large objects corresponding to a stiff EoS and would thus leave only a tiny
window of very massive stars in the mass-radius plane
\cite{Lavagetto:2006ew,Bejger:2006hn}
for a theory of compact star matter to fulfill all
above mentioned constraints.
\begin{figure} [ht]
\begin{tabular}{ll}
\includegraphics[angle=270,width=0.5\textwidth]{M-R-newest.ps}
\includegraphics[angle=270,width=0.5\textwidth]{dan_eos_newest.ps}&
\end{tabular}
\caption{Left panel: Mass and mass-radius constraints on compact star
configurations from recent observations compared to solutions of TOV equations
for nuclear EoS discussed in the text.
Right panel: Flow constraint from heavy-ion collisions
\cite{Danielewicz:2002pu} compared to the same set of nuclear EoS used in the
left panel. }
    \label{f:M-R}
\end{figure}

In the left panel of Fig.~\ref{f:M-R} we show some of these modern
observational constraints for masses and mass-radius relationships together
with solutions of the Tolman-Oppenheimer-Volkoff (TOV) equations for a set
of eight hadronic EoS classified in three groups:
(i) {\it relativistic mean-field (RMF) approaches}
{\it with non-linear (NL) self-interactions} of the $\sigma$
meson \cite{Gaitanos:2003zg}.
In NL$\rho$ the isovector part of the interaction is
described only by a $\rho$ meson, while the set NL$\rho\delta$ also includes a
scalar isovector meson $\delta$ that is usually neglected in RMF models
\cite{Liu02};
(ii) {\it RMF models with density dependent couplings and
masses}
are represented here by four different models from two classes, where in the
first one density dependent meson couplings are modeled so that a number of
properties of finite nuclei (binding energies, charge and diffraction radii,
surface thicknesses, neutron skin in ${}^{208}$Pb, spin-orbit splittings) can
be fitted \cite{Typel:2005ba}.
D${}^{3}$C has in addition a derivative coupling leading to momentum-dependent
nucleon self-energies and DD-F4 is modeled such that the flow constraint
\cite{Danielewicz:2002pu} from heavy-ion collisions is fulfilled.
The second class of these models is motivated by the Brown-Rho scaling
assumption \cite{Brown:1991kk} that not only the nucleon mass but also the
meson masses should decrease with increasing density.
In the KVR and KVOR models \cite{Kolomeitsev:2004ff}
these dependences were related to a
nonlinear scaling function of the  $\sigma$- meson field such that the EoS of
symmetric nuclear matter and pure neutron matter below four times the
saturation density coincide with those of the Urbana-Argonne group
\cite{APR}.
In this way the latter approach builds a bridge between the phenomenological
RMF models and (iii) {\it microscopic EoS} built on realistic
nucleon-nucleon forces.
Besides the variational approaches (APR \cite{APR}, WFF \cite{Wiringa:1988tp},
FPS \cite{Friedman:1981qw}) such ab-initio approaches to
nuclear matter are provided, e.g., by the relativistic
Dirac-\-Brueckner-\-Hartree-\-Fock (DBHF) \cite{DaFuFae04}
and the nonrelativistic Brueckner-\-Bethe-\-Goldstone
\cite{Baldo:1999rq} approaches.
Stiff EoS like D$^3$C, DD-F4, BBG and DBHF fulfill the demanding constraints
for a large radius and mass, while the softer ones like NL$\rho$ don't.
It is interesting to note that the flow constraint \cite{Danielewicz:2002pu}
shown in the right panel of  Fig.~\ref{f:M-R} sets limits to the tolerable
stiffness: it excludes the D$^3$C EoS and demonstrates that DD-F4, BBG and DBHF
become too stiff at high densities above $\sim 0.55$ fm$^{-3}$.
For a detailed discussion, see Ref. \cite{Klahn:2006ir}.

A key question asked in investigating the structure of matter at high
densities is how the quark substructure of hadrons manifests itself in the
EoS and whether the phase transition to quark matter can occur inside compact
stars.
In Ref. \cite{Ozel:2006km}, \"Ozel has debated that the new
constraints reported above would exclude quark matter in compact star
interiors reasoning that it would entail an intolerable
softening of the EoS. Alford et al. \cite{Alford:2006vz} have given a few
counter examples demonstrating that quark matter cannot be excluded.
In the following section we discuss a recently developed chiral quark model
\cite{Klahn:2006iw} which is in accord with the modern constraints,
see also \cite{Blaschke:2007ri}.

\section{Color superconducting quark matter: masquerade revisited}

We describe the thermodynamics of the deconfined quark matter phase
within a three-flavor quark model of Nambu--Jona-Lasinio (NJL) type, with a
mean-field thermodynamic potential given by
\begin{eqnarray}
\Omega_{MF}(T,\mu)
&=& \frac{1}{8 G_S}\left[\sum_{i=u,d,s}(m^*_i-m_i)^2
  - \frac{2}{\eta_V}(2\omega_0^2+\phi_0^2)
        +\frac{2}{\eta_D}\sum_{A=2,5,7}|\Delta_{AA}|^2\right]
        \nonumber \\
&&-\int\frac{d^3p}{(2\pi)^3}\sum_{a=1}^{18}
        \left[E_a+2T\ln\left(1+e^{-E_a/T}\right)\right]
        + \Omega_l - \Omega_0~.
\label{eos-blaschke-potential}
\end{eqnarray}
Here, $\Omega_l$ is the thermodynamic potential for electrons and muons,
and the divergent term $\Omega_0$ is subtracted in order to assure zero
pressure and energy density in vacuum ($T=\mu=0$). The quasiparticle
dispersion relations, $E_a(p)$, are the eigenvalues of the hermitean matrix
\begin{equation}
{\mathcal M} = \left[
\setlength\arraycolsep{-0.01cm}
\begin{array}{cc}
        -\vec{\gamma}\cdot\vec{p}- \hat{m}^*+\gamma^0\hat{\mu}^* &
        i \gamma_5 C \tau_A \lambda_A \Delta_{AA} \\
        i C \gamma_5 \tau_A \lambda_A\Delta_{AA}^* &
        -\vec{\gamma}^T\cdot\vec{p}+\hat{m}^*-\gamma^0\hat{\mu}^*
\end{array}
\right],
\label{eos-blaschke-eigmatrix}
\end{equation}
in color, flavor, Dirac, and Nambu-Gorkov space. Here, $\Delta_{AA}$ are the
diquark gaps. $\hat{m}^*$ is the diagonal renormalized mass matrix and
$\hat{\mu}^*$ the renormalized chemical potential matrix, $\hat{\mu}^*={\rm
  diag}_f
(\mu_u-G_S\eta_V\omega_0,\mu_d-G_S\eta_V\omega_0,\mu_s-G_S\eta_V\phi_0)$.
The gaps and the renormalized masses are determined by minimization of the
mean-field thermodynamic potential (\ref{eos-blaschke-potential}).
We have to obey constraints of charge neutrality which depend on the
application we consider.
In the (approximately) isospin symmetric situation of a heavy-ion collision,
the color charges are neutralized, while the electric charge in general is
non-zero.  For matter in $\beta$-equilibrium in compact stars, also the
global electric charge neutrality has to be fulfilled. For further details, see
\cite{Blaschke:2005uj,Ruster:2005jc,Abuki:2005ms,Warringa:2005jh}.

We consider $\eta_D$ as a free parameter of the quark matter model, to be
tuned with the present phenomenological constraints on the high-density EoS.
Similarly, the relation between the coupling in the scalar and vector meson
channels, $\eta_V$, is considered as a free parameter of the model.  The
remaining degrees of freedom are fixed according to the NJL model
parameterization in table I of \cite{Grigorian:2006qe}, where a fit to
low-energy phenomenological results has been made.

\begin{figure} [ht]
\begin{tabular}{ll}
\includegraphics[angle=270,width=0.5\textwidth]{M-R-Z_new.ps}
\includegraphics[angle=270,width=0.5\textwidth]{Flow_dbhf-njl.ps}
\end{tabular}
\caption{Same as Fig.~\ref{f:M-R} for hybrid EoS  with a low density hadronic
branch described by the DBHF approach and a high density quark matter branch
obtained from a three-flavor NJL model with color superconductivity (diquark
coupling $\eta_D$) and isoscalar vector meson coupling ($\eta_V$).}
    \label{f:P-n}
\end{figure}

As a unified description of quark-hadron matter, naturally
including a description of the phase transition, is not available yet,
although steps in this direction have been suggested within nonrelativistic
\cite{Ropke:1986qs} and relativistic \cite{Lawley:2006ps} models.
We apply here the so-called two-phase description, being aware of its
limitations. The nuclear matter phase is described within the
DBHF approach and the transition to the quark matter phase given above
is obtained by a Maxwell construction.
In the right panel of Fig.~\ref{f:P-n}
it can be seen that the necessary softening of the high density EoS in
accordance with the flow constraint is obtained for a vector coupling of
$\eta_V=0.5$ whereas an appropriate deconfinement density is obtained for a
strong diquark coupling in the range  $\eta_D=1.02 - 1.03$. The resulting phase
transition is weakly first order with an almost negligible density jump.
Applying this hybrid EoS with so defined free parameters under compact star
conditions a sequence of hybrid star configurations is obtained which fulfills
all modern constraints, see the left panel of Fig. 2.
In that figure we also indicate by a red dot the minimal mass $M_{DU}$ for
which the central density reaches a value allowing the fast direct Urca (DU)
cooling process in DBHF neutron star matter to occur, leading to problems
with cooling phenomenology \cite{Blaschke:2004vq,Blaschke:2006gd}.
Note that for a strong diquark coupling $\eta_D=1.03$, the critical density
for quark deconfinement is low enough to prevent the hadronic direct Urca (DU)
cooling problem by an early onset of quark matter. For the given hybrid
EoS, there is a long sequence of stable hybrid stars with two-flavor
superconducting (2SC) quark matter, before the occurrence of the strange quark
flavor and the simultaneous onset of the color-flavor-locking (CFL) phase
which renders the star gravitationally unstable
\cite{Baldo:2002ju,Klahn:2006iw}.
Comparing the hybrid star sequences with the purely hadronic DBHF ones one can
conclude that the former 'masquerade' themselves as neutron stars
\cite{Alford:2004pf} by having very similar mechanical properties.

Besides for mass-radius relationships, the masquerade effect can be discussed
also for the moment of inertia \cite{Klahn:2006iw}, which is becoming
accessible to measurements in relativistic binary systems like the double
pulsar J0737--3039, see \cite{Lattimer:2004nj,Bejger:2005jy}.
\begin{figure} [ht]
\includegraphics[width=7.5cm,height=14cm,angle=-90]{MG_MB_untrap.ps}
\caption{The effect of neutrino untrapping on hybrid star configurations:
a mass defect of $0.04~M_\odot$ occurs in the transition from neutrino
trapping case (solid red line, $Y_{le}=0.4$) to untrapping
(dashed black lines, $Y_{le}=0$) at fixed baryon mass $M_B$, independent
of the presence of a quark core for masses above the threshold indicated
by an asterix.
The blue rectangle corresponds to the constraint from the lighter
companion (star B) in the double pulsar PSR J0737-3039, see \cite{Pod05}.
}
\label{f:untrap}
\end{figure}

In these proceedings we add a new aspect to this masquerade discussion:the
energy release in the neutrino untrapping transition of a cooling protoneutron
star (PNS).
The process is pictured such that the PNS formed in a supernova collapse is
hot enough ($T_{\rm PNS}\sim 30 \dots 50$ MeV) to trap neutrinos since their
mean free path is shorter than the size of the star.
In this era the fast neutrino cooling cannot proceed from the volume but only
from the surface.
Until the neutrino opacity temperature ($T_{\rm opac} \sim 1$ MeV)
is reached, the neutrinos take part in the $\beta$-equilibration processes and
a finite lepton fraction is established. We will assume $Y_{le}=0.4$
corresponding to a neutrino chemical potential $\mu_\nu\sim 200$ MeV in the
PNS core.
For a discussion of hot hybrid PNS, see \cite{Nicotra:2006eg}.
In the untrapping transition at $T\sim T_{\rm opac}$, the neutrinos
decouple from the $\beta$-equilibrium so that $Y_{le}\to 0$ and the
equilibrium configuration gets readjusted with a new gravitational mass while
naturally, the baryon number is conserved in this process.
From Fig. \ref{f:untrap} one can estimate the mass defect (energy release) due
to neutrino untrapping to $\Delta M = 0.04~M_\odot$, almost independent of the
mass of the configuration and even of its structure: neutron stars and hybrid
stars result in the same energy release!
This is another aspect of the masquerade effect.
Note that the neutrino untrapped configurations fulfill the extended $M_G-M_B$
constraint \cite{Klahn:2006ir}, given by the rectangle in the right
panel of Fig.\ref{f:untrap}. This constraint is derived from the lighter
companion (star B) in the double pulsar J0737--3039 (B) \cite{Pod05}.

\section{Unmasking compact star interiors}
To unmask the neutron star interior
and to disentangle a quark matter core from a hadronic one,
might therefore require observables
based on transport properties, strongly modified from normal due to the color
superconductivity.
It has been suggested to base tests of the structure of matter at high
densities on analyses of the cooling behavior
\cite{Blaschke:2006gd,Popov:2004ey,Popov:2005xa}
or the stability of fastly rotating stars against r-modes
\cite{Madsen:1999ci,Drago:2007iy}.
It has turned out that for these phenomena the fine tuning of color
superconductivity in quark matter is an essential ingredient.
The result of these studies suggests that on the one hand unpaired quark matter
results in too fast cooling (same holds for pure two-flavor color
superconductivity (2SC) due to one unpaired (blue) color) and on the other,
complete pairing with large gaps (CFL phase, $\Delta_{\rm CFL}\sim 100$ MeV)
would result in too slow cooling \cite{Page:2000wt,Blaschke:2000dy},
r-mode instability \cite{Madsen:1999ci}
and gravitational collapse \cite{Baldo:2002ju,Klahn:2006iw}.
Viable alternatives for quark matter thus have to be color superconducting,
with all quark modes paired, but a few modes (at least one) should have a
small gap of $\Delta_X\sim 1$ MeV, desirable with a decreasing density
dependence, as suggested for the 2SC+X phase
\cite{Grigorian:2004jq,Grigorian:2006pu}.
Unfortunately, there is no microscopic model supporting the 2SC+X pairing
pattern yet. A suitable candidate, however, could be the recently suggested
single flavor pairing in the isotropic spin-1 color superconductivity phase
(iso-CSL) \cite{Aguilera:2005tg,Aguilera:2006cj}, see also
\cite{Marhauser:2006hy}.
It has been shown that the neutrino emissivity and bulk viscosity
of the iso-CSL phase fulfills constraints from cooling and r-mode stability
\cite{Blaschke:2007bv}.

\section{Conclusions}
In this contribution it is shown that modern quantum field
theoretical approaches to quark matter including color
superconductivity and a vector meanfield allow a microscopic description of
hybrid stars which fulfill the new, strong constraints
for high masses, large radii, sufficiently fast cooling evolution and
absence of r-mode instabilities.
The deconfinement transition in the resulting stiff hybrid equation of state
is weakly first order so that signals of it have to be expected due to
specific changes in transport properties governing the rotational and cooling
evolution caused by the color superconductivity of quark matter.
This conclusion holds analogously for the investigation of quark deconfinement
in future generations of nucleus-nucleus collision experiments
at low temperatures and high baryon densities,
such as CBM @ FAIR.
The extrapolation of the hybrid EoS constrained here from compact star physics
results in almost crossover behavior at the deconfinement transition
 \cite{Grigorian:2006pu}, so that promising observable effects for its
diagnostics shall be based on (precursor) effects of color superconductivity.

\section*{Acknowledgements}
We thank the organizers of the EXOCT07 conference for providing a perfect
environment for discussios of  exotic matter and compact star research.
D.B. is supported by the Polish Ministry of Science and Higher
Education, T.K. is grateful for partial support from GSI Darmstadt and the
Department of Energy, Office of Nuclear Physics, contract no.\
DE-AC02-06CH11357.
F.S. acknowledges support from the Swedish Graduate
School of Space Technology and the Royal Swedish Academy of Sciences.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction} \label{sec:intro}

 In this paper we consider the following two-dimensional elliptic optimal control
problem  \cite{rees2010optimal}: Find the state $u^{*} \in H^{1}(\Omega)$ and the control $f^{*}\in L^{2}(\Omega)$ such that
\begin{equation}\label{eq:Control-problem}
  \mathcal{J}(u^{*},f^{*}) = \displaystyle\min_{(u,f)\in H^{1}(\Omega)\times L^{2}(\Omega)}\mathcal{J}(u,f),
\end{equation}
subject to the state equation
\begin{equation}\label{eq:sub-state-equ}
\left\{
  \begin{aligned}
  -\triangle u &=&f,\,\,\,\, {\rm in}\,\,\Omega,\\
  u&=&g,\,\,\, {\rm at}\,\,\, \partial \Omega,\\
   \end{aligned}
               \right. \\
\end{equation}
with cost functional
\begin{equation}\label{eq:cost-fun}
  \mathcal{J}(u,f) = \frac{1}{2}\|u-\hat{u}\|^2_{L^2(\Omega)}+\beta\|f\|^2_{L^2(\Omega)},
\end{equation}
where $\hat{u}$ is the desired state and $\beta>0$ is the weight of the cost of the control
(or simply a regularization parameter).

The discretization of the optimization problem \eqref{eq:Control-problem} leads to a large scale saddle-point system.  Since  multigrid methods offer
the possibility of solving problems with $N$ unknowns using $O(N)$ work and storage, and lead to substantial improvements in computational efficiency over direct methods, they have gained growing interest in the area of optimal control problems, see,  for example \cite{borzi2003multigrid,borzi2007high,borzi2009multigrid, engel2011multigrid, schoberl2011robust, simon2009schwarz, takacs2011convergence}.

The choice of multigrid smoother is crucial to construct efficient multigrid algorithms.  Local Fourier analysis (LFA) is a useful tool to  help choose multigrid components, such as relaxation schemes and grid-transfer operators, and can quantitatively predict the two-grid convergence factor.  The main goal of this work is to construct and analyze  multigrid methods  for the discrete optimal control problems.  Note that there are a few studies of LFA applied to the optimal control problems.  In \cite{kepler2009fourier},  the smoothing factor for a one-dimensional optimal control problem was proposed by LFA for a special smoother. \cite{borzi2009multigrid} reviewed some recent efforts and outlined recent developments in the field of multigrid methods for partial differential equation optimization, where LFA was applied to study collective Gauss–Seidel methods. The finite difference multigrid solution of an optimal control problem associated with an elliptic equation was  considered in \cite{borzi2002accuracy}, and LFA was used to estimate the convergence factor for the collective Gauss–Seidel relaxation.

 In this work,  we consider  finite element  monolithic multigrid solution of the  elliptic distributed control problem  \eqref{eq:Control-problem}, that is, employing multigrid methods directly to the discrete system in coupled form.   In contrast to collective smoothers, we consider a well-known Braess-Sarazin-type relaxation (BSR) scheme, which was original  designed for the Stokes equations \cite{braess1997efficient}.   In our recently work \cite{CH2021addVanka}, we have proved that the mass matrix  obtained from bilinear finite elements is a good approximation to the inverse of the scalar Laplacian discretized by a five-point finite difference method, where the optimal smoother factor of $\frac{1}{3}$ is derived for the Laplacian. Furthermore, we extended this mass-based approximation to a Braess-Sarazin-type relaxation scheme for the Stokes equations \cite{YH2021massStokes}, where we obtained the optimal smoother factor of $\frac{1}{3}$  for Stokes. This motivates us to consider the stiffness matrix derived from the finite-difference discrerization of  the scalar Laplacian to approximate the mass matrix arsing from the finite element discretizations of optimal control problems.  In the resulting BSR relaxation, we can directly apply the stiffness matrix and there is no need to invert a matrix. LFA is applied  to examine this stiffness-based BSR applied to the optimal control problems considered here. We derive  optimal smoothing factor of $\frac{1}{3}$ for the stiffness-based BSR for the optimal control problems. Considering practical use, an inexact version of stiffness-based Braess-Sarazin relaxation is developed, where a  3 $V(2,2)$-cycles  multigrid  with weighted Jacobi relaxation is used to approximate the solution of Schur complement system appeared in the exact BSR process.  This inexact version preserves the optimal smoothing factor of $\frac{1}{3}$ obtained from exact BSR. The LFA  convergence estimates agree very well with results of numerical experiments and   are independent of the grid size and of the value of the control parameter.

This paper is organized as follows.  In Section \ref{sec: Discretization-relaxation}, we propose a  stiffness-based Braess-Sarazin relaxation scheme for the  elliptic distributed control problem.
In Section \ref{sec:smoothing-analysis}, we apply LFA to the proposed stiffness-based Braess-Sarazin relaxation, and we obtain a highly satisfactory convergence factor. Numerical results are reported in Section \ref{sec:Numer} to confirm our theoretical analysis. Finally, we draw conclusions in Section \ref{sec:concl}.

\section{Discretization and relaxation}\label{sec: Discretization-relaxation}

In our work, we  discretize  problem \eqref{eq:Control-problem}  using bilinear quadrilateral $Q_1$ finite elements. The discretize formulations of \eqref{eq:sub-state-equ} and \eqref{eq:cost-fun} are given by
\begin{eqnarray}
&&\min_{u_h,f_h}\frac{1}{2}\|u_h-\hat{u}\|^2_{L^2(\Omega)}+\beta\|f_h\|^2_{L^2(\Omega)},\label{eq:weak-form1}\\
&&{\rm such\,\, that} \int \nabla u_h v_h = \int v_h f_h, \,\,\forall v_h\in V_0^h, \label{eq:weak-form2}
\end{eqnarray}
where $V_0^h$ is the $Q_1$ finite element space. With a Lagrange multiplier ($\tau$) approach to  \eqref{eq:weak-form1} and  \eqref{eq:weak-form2}, it typically leads to a linear system of the form  \cite{rees2010optimal}
 \begin{equation}\label{eq:3-block-system}
\begin{pmatrix}
2 \beta M   &     0     & -M^T\\
0              & M       &K^T \\
 -M           &K       & 0
 \end{pmatrix}
\begin{pmatrix} f_h \\ u_h \\ \tau_h \end{pmatrix}
 =b_h,
 \end{equation}

 where $K$ is  the stiffness matrix (the discrete Laplacian) and  $M$ is the discrete mass matrix. For more details, we refer to \cite{rees2010optimal}.

We can rewrite  system  \eqref{eq:3-block-system} as a saddle-point system, that is,

 \begin{equation}\label{eq:saddle-structure}
       \mathcal{L}_{h}  z=\begin{pmatrix}
      A & B^{T}\\
     B & 0\\
    \end{pmatrix}
        \begin{pmatrix} x \\ y
        \end{pmatrix}
  =b_h,
 \end{equation}

where
\begin{equation*}
A=\begin{pmatrix}
2 \beta M   &     0 \\
0              & M
\end{pmatrix},\quad
B= \begin{pmatrix}
-M  &     K
\end{pmatrix}.
\end{equation*}
Here, we consider multigrid methods for solving \eqref{eq:saddle-structure}. One important process in multgrid is the smoothing, where we update a current approximation $z_k$ via
\begin{equation}\label{eq:general relaxation-form}
\hat{z}_{k} = z_k-\omega \mathcal{K}^{-1}_h (b_h-  \mathcal{L}_{h}  z_k)=z_k-\omega \delta z,
\end{equation}
where $\mathcal{K}^{-1}_h$ is an approximation to  $\mathcal{L}_{h} $, and $\delta z =\mathcal{K}^{-1}_h (b_h  - \mathcal{L}_{h} z_k)$.   Then, the error-propagation for the relaxation scheme \eqref{eq:general relaxation-form} is
\begin{equation}\label{eq:Relaxation-error-S}
\mathcal{S}_h=I-\omega   \mathcal{K}^{-1}_h \mathcal{L}_h.
\end{equation}
The choices of $\mathcal{K}_h$ and $\omega$  play a crucial role in determining the convergence speed of multigrid methods.   For
simplicity, throughout the rest of this paper, we drop the subscript $h$, except when necessary for clarity.

In this work, we consider  Braess-Sarazin relaxation \cite{braess1997efficient}, that is,   $\mathcal{K}$ has the form of
\begin{equation}\label{eq:Precondtion}
   \mathcal{K}=  \begin{pmatrix}
    \alpha C & B^{T}\\
     B & 0\\
    \end{pmatrix},
\end{equation}
where $C$ is an approximation of $A$, whose inverse is easy to apply, and $\alpha>0$.

With $\mathcal{K}$ defined in \eqref{eq:Precondtion}, we compute the update $\delta z=(\delta x, \delta y)$ in \eqref{eq:general relaxation-form} by  the following two stages
\begin{eqnarray}
  (BC^{-1}B^{T})\delta y&=&BC^{-1}r_{x}-\alpha r_{y}, \label{eq:solution-of-precondtion-s1}\\
  \delta x&=&\frac{1}{\alpha}C^{-1}(r_{x}-B^{T}\delta y),\label{eq:solution-of-precondtion-s2}
\end{eqnarray}
where $(r_{x},r_{y})=b_h-\mathcal{L}z_k$.

Note that, in \eqref{eq:solution-of-precondtion-s1} and   \eqref{eq:solution-of-precondtion-s2}, we can directly consider an approximation to  $C^{-1}$.   In our recently work \cite{CH2021addVanka}, we have shown that the mass matrix $M$ derived from bilinear finite elements is a good approximation to the inverse of the scalar Laplacian discretized by the five-point finite difference method. This inspires us to  consider the  discrete scalar Laplacian approximating  the inverse of the mass matrix $M$ in  the diagonal part of  the coefficient matrix defined in \eqref{eq:3-block-system}. Let $A_{fd}$ be the stiffness matrix obtained from the five-point finite difference discretization for the scalar Laplacian. In \eqref{eq:Precondtion}, we consider
\begin{equation}\label{eq:C-inverse-Q}
C^{-1}=C^{-1}_f =
\begin{pmatrix}
\frac{1}{2\beta}A_{fd} &  0\\
0 &  A_{fd}
\end{pmatrix},
\end{equation}
and denote the corresponding  $\mathcal{K}$ as
\begin{equation}\label{eq:Precondtion-K1}
   \mathcal{K}_f=  \begin{pmatrix}
    \alpha C_f& B^{T}\\
     B & 0\\
    \end{pmatrix}.
\end{equation}
To choose suitable values of  the parameters $\omega$ and $\alpha$, we consider LFA to examine the relaxation scheme with $\mathcal{K}=\mathcal{K}_f$ in the following section.

\section{Smoothing analysis}\label{sec:smoothing-analysis}
LFA \cite{MR1807961, wienands2004practical} has been widely used to analyze  actual multigrid performance for different types of problems, by investigating  the spectral radius of {\it symbol} of the underlying operator.  Specifically, the LFA smoothing factor  with simplifying  assumptions on boundary conditions offers a  sharp convergence estimate of multigrid methods. Thus, in this work, we employ LFA to examine the LFA smoothing factor for  BSR,  and help choose relaxation parameters.  Next, we first give a brief introduction to LFA.

We consider two-dimensional infinite uniform grids $\mathbf{G}_{h}$
\begin{equation*}
  \mathbf{G}_{h}=\left\{\boldsymbol{x}:=(x_1,x_2)=\boldsymbol{k}h=(k_{1},k_{2})h, k\in \mathbb{Z}\right\},
\end{equation*}
and Fourier modes $\varphi(\boldsymbol{\theta},\boldsymbol{x})=e^{\imath \boldsymbol{\theta}\boldsymbol{x}/h}$ on $\mathbf{G}_{h}$, where $\boldsymbol{\theta}=(\theta_1,\theta_2)$ and $\imath^2=-1$.  Since $\varphi(\boldsymbol{\theta},\boldsymbol{x})$ is periodic in $\boldsymbol{\theta}$ with period $2\pi$, we consider $\theta_i\in\big(-\frac{\pi}{2},\frac{3\pi}{2}\big]$. The coarse grid $  \mathbf{G}_{2h}$ is defined similarly.

Let $L_h$ be a Toeplitz operator acting on $\mathbf{G}_{h}$ as follows \cite{MR1807961},
\begin{eqnarray*}
  L_{h}w_{h}(\boldsymbol{x})&=&\sum_{\boldsymbol{\kappa}\in\boldsymbol{V}}s_{\boldsymbol{\kappa}}w_{h}(\boldsymbol{x}+\boldsymbol{\kappa}h),
  \end{eqnarray*}
with constant coefficients $s_{\boldsymbol{\kappa}}\in \mathbb{R}$, where $w_{h}(\boldsymbol{x})$ is a function in $l^{2} (\mathbf{G}_{h})$. Here, $\boldsymbol{V}$ is a finite index set.

\begin{definition}\label{formulation-symbol}
  We call $\widetilde{L}_{h}(\boldsymbol{\theta})=\displaystyle\sum_{\boldsymbol{\kappa}\in\boldsymbol{V}}s_{\boldsymbol{\kappa}}e^{\iota \boldsymbol{\theta}\cdot \boldsymbol{\kappa}}$ the symbol of $L_{h}$.
\end{definition}

We consider multigrid methods for finite-element discretizations with standard geometric grid coarsening, that is, we construct a sequence of coarse grids by doubling the mesh size in each spatial direction.  High and low frequencies for standard coarsening are given by
\begin{equation*}
  \boldsymbol{\theta}\in T^{{\rm low}} =\left[-\frac{\pi}{2},\frac{\pi}{2}\right)^{2}, \, \boldsymbol{\theta}\in T^{{\rm high}} =\displaystyle \left[-\frac{\pi}{2},\frac{3\pi}{2}\right)^{2} \bigg\backslash \left[-\frac{\pi}{2},\frac{\pi}{2}\right)^{2}.
\end{equation*}

\begin{definition} The LFA smoothing factor for  the error-propagation for the relaxation scheme $\mathcal{S}$, see \eqref{eq:Relaxation-error-S}, is defined as
\begin{equation}\label{eq:LFA-mu-form}
  \mu_{{\rm loc}}(\boldsymbol{p})=\max_{\boldsymbol{\theta}\in T^{{
  \rm high}}}\big\{\big|\lambda(\widetilde{\mathcal{S}}(\boldsymbol{\theta}))\big| \,\,\big\},
\end{equation}
where  $\boldsymbol{p}$ is algorithmic parameters,  and $\lambda\big(\widetilde{\mathcal{S}}(\boldsymbol{\theta})\big)$ denotes the  eigenvalue of symbol $\widetilde{\mathcal{S}}(\boldsymbol{\theta})$.
\end{definition}
Note that for the $Q_1$  elements considered in this work,  $\widetilde{\mathcal{S}}(\boldsymbol{\theta})$ is a $3\times 3$ matrix. Since $\mu_{{\rm loc}}(\boldsymbol{p})$ is a function of   $\boldsymbol{p}$, we can minimize  $\mu_{{\rm loc}}(\boldsymbol{p})$ to obtain a fast convergence speed.  Thus, we define the optimal smoothing factor as follows.
\begin{definition}
The optimal LFA  smoothing factor  for  the error-propagation for the relaxation scheme $\mathcal{S}$ is defined as
  \begin{equation*}
    \mu_{{\rm opt}}=\min_{\boldsymbol{p} \in \Upsilon}{\mu_{{\rm loc}}}(\boldsymbol{p}),
  \end{equation*}
  where $\Upsilon$ is the set of allowable parameters.
\end{definition}

Another important LFA factor is the LFA two-grid convergence factor.  In general, the two-grid error-propagation can be expressed as
\begin{equation*}
E_h = \mathcal{S}^{\nu_2} (I-P (\mathcal{L}_{2h})^{-1} R \mathcal{L}_h)\mathcal{S}^{\nu_1},
\end{equation*}
where $\mathcal{L}_{2h}$ is the coarse grid operator,  $R$ and $P$ are the restriction and interpolation operators, respectively.  The integers $\nu_1$ and $\nu_2$ are pre- and post- smoothing steps, respectively.

\begin{definition}
The   LFA  two-grid convergence factor for $E_h$ is defined as
  \begin{equation}\label{eq:LFA-rho-form}
    \rho=\max_{\boldsymbol{\theta}\in T^{{
  \rm low}}} \left\{\rho(\widetilde{E}_h(\boldsymbol{\theta}))\right\} ,
  \end{equation}
  where $\widetilde{E}_h(\boldsymbol{\theta})$ is  the two-grid LFA symbol of  $E_h$, and $\rho(\widetilde{E}_h)$ denotes the spectral radius of matrix $\widetilde{E}_h$.
\end{definition}
Note that in our case,  $\widetilde{E}_h(\boldsymbol{\theta})$  is a $12\times 12$ matrix. For more details on how to compute two-grid symbol, we refer to  \cite{MR1807961}. Since it is easy to compute the LFA smoothing factor, \eqref{eq:LFA-mu-form}, compared to the LFA two-grid convergence factor, \eqref{eq:LFA-rho-form}, we focus on the analysis of LFA smoothing factor in  the following.

 \subsection{Fourier representation of mass and stiffness operators}

In one dimension (1D), the  stiffness and mass stencils of $Q_1$ elements are
\begin{equation}\label{eq:1D-Q1Q1-stencil}
K_1=\frac{1}{h}\begin{bmatrix}  -1 & 2 &-1  \end{bmatrix} \quad \text{and} \quad
M_1=\frac{h}{6}\begin{bmatrix}  1 & 4 &1  \end{bmatrix},
\end{equation}
respectively.

Then, the stiffness stencil for bilinear discretization in two dimensions (2D)  is given by
\begin{eqnarray*}
  K &=& K_1\otimes M_1+M_1\otimes K_1 \nonumber\\
   &=& \frac{1}{3}\begin{bmatrix}
   -1 & -1 &-1\\
   -1 & 8 &-1\\
   -1 & -1 &-1
   \end{bmatrix},
\end{eqnarray*}
and the mass stencil for bilinear discretization  in 2D is given by
\begin{eqnarray}
  M &=& M_1\otimes M_1 \nonumber\\
   &=& \frac{h^2}{36}\begin{bmatrix}
   1 &  4    & 1\\
   4 &  16   & 4\\
   1 &  4    & 1
   \end{bmatrix}.\label{eq:mass-2D-stencil}
\end{eqnarray}
From  \eqref{eq:1D-Q1Q1-stencil}, the symbols of the stiffness and mass stencils in 1D are
\begin{equation}\label{eq:1D-stencil-symbol}
\widetilde{K}_1(\theta)=\frac{2}{h}(1-\cos\theta),\,\,\,
\widetilde{M}_1(\theta)=\frac{h}{3}(2+\cos\theta).
\end{equation}
Based on \eqref{eq:1D-stencil-symbol}, the symbol of $K$ can be written as
\begin{eqnarray*}
  \widetilde{K}(\theta_1,\theta_2) &=&\widetilde{K}_1(\theta_2)\widetilde{M}_1(\theta_1)+
  \widetilde{M}_1(\theta_2) \widetilde{K}_1(\theta_1) \nonumber \\
  &=&\frac{2}{3}(4-\cos\theta_1-\cos\theta_2-2\cos\theta_1\cos\theta_2),
\end{eqnarray*}
and the symbol of $M$ is
\begin{eqnarray*}
  \widetilde{M}(\theta_1,\theta_2) &=&\widetilde{M}_1(\theta_2) \widetilde{M}_1(\theta_1)\nonumber \\
  &=&\frac{h^2}{9}(4+2\cos\theta_1+2\cos\theta_2+\cos\theta_1\cos\theta_2).
\end{eqnarray*}
We introduce  the stencil of the Laplacian $-\triangle_{h}$ discretized by the standard five-point scheme, given by
\begin{equation*}
 A_{fd}=  -\triangle_{h} =\frac{1}{h^2}\begin{bmatrix}
       & -1 &  \\
      -1 & 4 & -1 \\
        & -1 &
    \end{bmatrix}.
 \end{equation*}

Then,  the symbol of $A_{fd}$ is
\begin{equation}\label{eq:scalar-symbol-Laplace}
\widetilde{A}_{fd} =\frac{4-2\cos\theta_{1}-2\cos\theta_2}{h^2}.
\end{equation}
For simplicity, let $a= \widetilde{M}$ and $\hat{a}=(\widetilde{A}_{fd})^{-1}$.
Using \eqref{eq:scalar-symbol-Laplace}   gives
\begin{equation} \label{eq:AM-MA-symbol}
\widetilde{A}_{fd} \widetilde{M} =\widetilde{M} \widetilde{A}_{fd}= \frac{a}{\hat{a}}.
\end{equation}

\subsection{LFA for Braess-Sarazin relaxation scheme}\label{subsec:BSR-analytical-result}

Now, we examine  Braess-Sarazin-type algorithms, which  were originally developed as a relaxation scheme for the Stokes equations \cite{braess1997efficient}.   In practice,  \eqref{eq:solution-of-precondtion-s1} is not solved exactly. Thus, in the following, we first consider  exact solve for  \eqref{eq:solution-of-precondtion-s1}, and then inexact one.

 {\bf Exact Braess-Sarazin relaxation}: In \eqref{eq:Precondtion},   we consider $C^{-1}$  defined in \eqref{eq:C-inverse-Q} and $\alpha=1$. We first derive the corresponding optimal smoothing factor for the exact Braess-Sarazin relaxation.   We  refer to the relaxation scheme as  stiffness-based BSR.
Let $b= \widetilde{K}(\theta_1,\theta_2) $.
The symbol of operator $\mathcal{L}$ defined in \eqref{eq:saddle-structure} is given by
\begin{equation*}
  \widetilde{\mathcal{L}}(\theta_1,\theta_2) = \begin{pmatrix}
       2\beta a&       0 &          -a  \\
                0 &      a &           b \\
              -a  &       b&           0
    \end{pmatrix},
\end{equation*}
and
\begin{equation*}
  \widetilde{\mathcal{K}}_f(\theta_1,\theta_2) =\begin{pmatrix}
       2\beta \hat{a}&       0 &          -a  \\
                0 &       \hat{a} &           b \\
              -a  &       b&           0
    \end{pmatrix}.
\end{equation*}
Then,
\begin{equation*}
  \widetilde{\mathcal{L}}-  \lambda \widetilde{\mathcal{K}}_f
   = \begin{pmatrix}
   2\beta (a-\lambda\hat{a})&              0                     &          -a(1-\lambda)  \\
                0 &                      a-\lambda \hat{a}         &           b(1-\lambda) \\
 -a(1-\lambda)  &                   b(1-\lambda)                            &           0
    \end{pmatrix}.
\end{equation*}
The determinant of $\widetilde{\mathcal{L}}-  \lambda \widetilde{\mathcal{K}}_f$ is

\begin{align*}
 |\widetilde{\mathcal{L}}-  \lambda \widetilde{\mathcal{K}}_f | &= -2\beta (a-\lambda\hat{a}) b^2(1-\lambda)^2-a^2(1-\lambda)^2(a-\lambda \hat{a})\\
  & =(1-\lambda)^2 (a-\lambda \hat{a})(-2\beta b^2-a^2).
\end{align*}
It follows that the eigenvalues of   $\widetilde{\mathcal{K}}^{-1}_f\widetilde{\mathcal{L}}$ are $1,1$ and $\frac{a}{\hat{a}}$.

We restate the results for  solving the Laplace problem using  mass matrix approximation \cite{CH2021addVanka} in the following.

\begin{lemma} \cite{CH2021addVanka}\label{lem:ratio-a-hat-a}
For $\boldsymbol{\theta} \in T^{high}$, $\frac{a}{\hat{a}}(\boldsymbol{\theta}) \in [8/9, 16/9]$.
\end{lemma}
\begin{lemma} \cite{CH2021addVanka}\label{lemma:mass-smoothing-factor}
If we consider the mass matrix $M$ to approximate the inverse of the scalar  Laplacian, $A_{fd}$, and  $\mathcal{S}=I-\omega MA_{fd}$, then  the optimal smoothing factor for $\mathcal{S}$  is
\begin{equation*}
  \mu_{\rm opt}= \min_{\omega} \max_{\boldsymbol{\theta}\in T^{{
  \rm high}}} \left|1- \omega \widetilde{ M} \widetilde{A}_{fd}\right|=\frac{1}{3},
\end{equation*}
where the minimum is uniquely achieved at $\omega=\omega_{\rm opt} = \frac{3}{4}$.
\end{lemma}

Since $\widetilde{ M}\widetilde{A}_{fd}=\widetilde{A}_{fd}\widetilde{ M}=\frac{a}{\hat{a}}$, see \eqref{eq:AM-MA-symbol}, we have the following result.
\begin{corollary}\label{corollary:laplace-approximate-mass}
If we consider the discrete Laplacian $A_{fd}$ to approximate the inverse of the mass matrix, $M$, and  $\mathcal{S}_{fd}=I-\omega A_{fd} M$, then  the optimal smoothing factor for $\mathcal{S}_{fd}$  is
\begin{equation*}
  \mu_{\rm opt}= \min_{\omega} \max_{\boldsymbol{\theta}\in T^{{
  \rm high}}}  \left|1- \omega \widetilde{A}_{fd}\widetilde{ M}\right| =\frac{1}{3},
\end{equation*}
provided that $\omega=\omega_{\rm opt} = \frac{3}{4}$.
\end{corollary}
Since  $\widetilde{\mathcal{K}}^{-1}_f\widetilde{\mathcal{L}}$ contains eigenvalue $\frac{a}{\hat{a}}$, which is an eigenvalue of $\widetilde{A}_{fd}\widetilde{ M}$, Corollary \ref{corollary:laplace-approximate-mass}  indicates that  $\frac{1}{3}$ is a lower bound on the optimal smoothing factor for the  stiffness-based  BSR.
\begin{theorem}\label{thm:EBSR-smoothing-theorem}
The optimal smoothing factor for  the stiffness-based BSR for the optimal control system \eqref{eq:saddle-structure} is
  \begin{equation*}
   \mu_{{\rm opt},B}=\displaystyle \min_{\omega_{B}}\max_{ \boldsymbol{\theta}\in T^{{\rm high}}} \left|\lambda( \widetilde{\mathcal{S}}_{B}( \omega_{B},\boldsymbol{\theta}))\right|=\frac{1}{3},
\end{equation*}
where the minimum is attained provided that $\omega_{B}=\omega_{B,\rm opt}=\frac{3}{4}$.
\end{theorem}
\begin{proof}
 Since the eigenvalues of $\widetilde{\mathcal{S}}_{B}=I-\omega_B\widetilde{\mathcal{K}}^{-1}_f\widetilde{\mathcal{L}}$ are $1-\omega_B, 1-\omega_B$ and $1-\omega_B\frac{a}{\hat{a}}$ and Lemma \ref{lem:ratio-a-hat-a}, the smoothing factor for $\mathcal{\widetilde{S}}_{B}$ is
 \begin{equation*}
\mu_{{\rm loc}} = \max \left\{|1-\omega_B|, \left|1-\frac{16}{9}\omega_B\right|, \left|1-\frac{8}{9}\omega_B\right|\right \},
\end{equation*}
From Corollary \ref{corollary:laplace-approximate-mass}, we know that
 \begin{equation*}
\min_{\omega_B}\max \left\{\left |1-\frac{16}{9}\omega_B\right|, \left|1-\frac{8}{9}\omega_B\right|\right \}=\frac{1}{3},
\end{equation*}
provided that $\omega_B=\omega_{B,\rm opt}=\frac{3}{4}$.  Furthermore,    $|1-\omega_{B,\rm opt}|=\frac{1}{4}<\frac{1}{3}$.   Thus,  $\min_{\omega_B}\mu_{{\rm loc}} =\frac{1}{3}$.
\end{proof}
 \begin{remark}
 We point out that one can consider $\alpha$ as a variable in \eqref{eq:Precondtion-K1} rather than $\alpha=1$. However,  this does not change the optimal smoothing factor  for exact BSR obtained in Theorem \ref{thm:EBSR-smoothing-theorem}, and only changes the choice of optimal parameter $\omega_B$. In contrast to exact BSR,  varying $\alpha$  is needed for inexact version of BSR, which will be discussed in Section \ref{sec:Numer}.
 \end{remark}

 {\bf Inexact Braess-Sarazin  relaxation}:  Exact  Braess-Sarazin algorithm requires an exact inversion of the Schur complement \eqref{eq:solution-of-precondtion-s1}, which is very expensive.  In the literature, it has been shown that   inexact  solve  using a few sweeps of weighted Jacobi iteration for the Schur complement system is enough to maintain the same convergence of exact BSR for the Stokes equations, see, for example \cite{YH2021massStokes}.   In \eqref{eq:solution-of-precondtion-s1}, the Schur complement matrix is
 \begin{align*}
 BC^{-1} B^T &= \begin{pmatrix}
-M & K
\end{pmatrix}
\begin{pmatrix}
\frac{1}{2\beta}A_{fd} &  0\\
0 &  A_{fd}
\end{pmatrix}
 \begin{pmatrix}
-M^T \\   K^T
\end{pmatrix} \\
& = \frac{1}{2\beta} MA_{fd}M^T+KA_{fd}K^T.
\end{align*}
We explore  inexact stiffness-based BSR, using either a few sweeps of weighted Jacobi iteration to approximate the solution of \eqref{eq:solution-of-precondtion-s1} or a $V$-cycles multigrid on the Schur complement system. We are wondering whether the inexact version is possible to achieve the same smoothing factor of $\frac{1}{3}$. Here, we do not  seek theoretical analysis. Instead, we numerically examine  the performance of  inexact BSR in Section \ref{sec:Numer}.

Note that another choice for $\mathcal{K}$ is
\begin{equation}\label{eq:KJ-form}
\mathcal{K}=\mathcal{K}_D=\begin{pmatrix}
C_D & B^T\\
B  &  0
\end{pmatrix},
\end{equation}
where
\begin{equation*}
C_D=\begin{pmatrix}
2\beta {\rm diag}(M)  & 0\\
0  & {\rm diag}(M)
\end{pmatrix}.
\end{equation*}
From \eqref{eq:mass-2D-stencil}, it can be seen that the symbol of ${\rm diag}(M)$ is $\frac{4h^2}{9}=:a_D$.
Here, we are curious about the optimal smoothing factor for the smoother  \eqref{eq:KJ-form}. As a comparison, we derive the optimal smoothing factor for this choice.

\begin{theorem}\label{thm:Jacobi-BSR-smoothing}
If we consider $\mathcal{K}_D$ defined in \eqref{eq:KJ-form} and $\mathcal{S}_D=I-\omega_D \mathcal{K}^{-1}_D\mathcal{L}$,  them  the optimal smoothing factor for $\mathcal{S}_D$ is
 \begin{equation*}
   \mu_{{\rm opt},D}=\min_{\omega_{D}}\max_{ \boldsymbol{\theta}\in T^{{\rm high}}} \left|\lambda(\widetilde{\mathcal{S}}_{D}( \omega_{D},\boldsymbol{\theta}))\right|=\frac{5}{7},
\end{equation*}
provided that $\omega_D =\frac{8}{7}$.
\end{theorem}

\begin{proof}
It can easily be shown that the eigenvalues of $\widetilde{\mathcal{S}}_D$ are $1,1$ and $\frac{a}{a_D}$. Now, we consider the optimal smoothing factor for the modes, $\frac{a}{a_D}$. Note that

\begin{equation*}
\frac{a}{a_D} = \frac{1}{4}(4+2\cos\theta_1+2\cos\theta_2+\cos\theta_1\cos\theta_2).
\end{equation*}
For $\boldsymbol{\theta}\in T^{\rm high}$, it can easily be shown that the extreme values of $\frac{a}{a_D}$ are $\frac{1}{4}$ and $\frac{3}{2}$. Thus,

 \begin{equation*}
\mu_{\rm opt}=\min_{\omega_D}\max_{ \boldsymbol{\theta}\in T^{{\rm high}}} \left\{  |1- \omega_D|,\left |1-\frac{1 }{4}\omega_D\right|, \left|1-\frac{3}{2}\omega_D\right|\right \}=\frac{5}{7},
\end{equation*}
provided that $\omega_D=\frac{2}{1/4+3/2}=\frac{8}{7}$.
\end{proof}

From Theorems \ref{thm:EBSR-smoothing-theorem} and   \ref{thm:Jacobi-BSR-smoothing}, we see that the stiffness-based BSR offers a highly satisfactory convergence speed.  Thus, in our numerical tests, we only consider the stiffness-based BSR.

\section{Numerical experiments}\label{sec:Numer}
In this section, we  first present the LFA two-grid convergence to compare with the optimal LFA smoothing factor derived in Subsection \ref{subsec:BSR-analytical-result} for the exact stiffness-based BSR.  Second,  we display actual multigrid performance to validate our analytical result for exact BSR. To consider practical use, we  develop an inexact version of BSR, where a $V$-cycle multigrid is applied to the Schur complement system to maintain the performance observed for exact BSR.  To demonstrate the robustness of our multigrid methods with respect to the regularization parameters, we consider different values of $\beta$ in our tests.

\subsection{LFA predictions}

In this subsection, we report the LFA predicted two-grid convergence factor, $\rho_h$, for different values of regularization parameter $\beta$ and grid size $h$ for exact stiffness-based BSR.   For the prolongation of correction,  we consider the standard bilinear  interpolation and apply the corresponding adjoint operator for the restriction.  Galerkin operator is used to define the coarse grid operator. The number of smoothing steps $\nu$ is   $\nu = \nu_1+\nu_2$. We use $32 \times 32$ evenly distributed Fourier frequencies in the Fourier domain $[-\frac{\pi}{2}+\tau, \frac{\pi}{2}-\tau]^2$ with $\tau=\frac{\pi}{64}$ to compute LFA two-grid convergence factor  $\rho_h$ in \eqref{eq:LFA-rho-form}.

In Table \ref{tab:LFA-beta-rho}, we show the LFA predicted two-grid convergence factor as a function of $\beta$ and $h$ with $\nu=1$. We see that the LFA  two-grid convergence factor is the same as the optimal smoothing factor derived from our theoretical analysis. Moreover,  the convergence factor is robust in  regularization parameter $\beta$ and is $h$-independent.

 \begin{table}[H]
 \caption{LFA predicted two-grid convergence factor, $\rho_h$,  as a  function of $\beta$ and $h$ with $\nu=1$ for exact stiffness-based BSR.}
\centering
\begin{tabular}{ l    c c c c}
\hline
$h$       & $\beta=10^{-2}$  &$\beta=10^{-4}$    &$\beta=10^{-6}$   & $\beta=10^{-8}$  \\
\hline
 $1/64$                    &0.333        &0.333       &0.333      &0.333   \\
 $1/128$                  &0.333        &0.333       &0.333      &0.333  \\
 $1/256$                   &0.333        &0.333       &0.333      &0.333    \\
\hline
\end{tabular}\label{tab:LFA-beta-rho}
\end{table}
In Table \ref{tab:rho-nu}, we report the LFA predicted two-grid convergence factor as a function of the number of smoothing steps $\nu$  and $h$ with $\beta=10^{-2}$. We see that $\rho_h=\mu^{\nu}_{\rm opt}$ for different $h$.  For other $\beta$, we see same behaviour,  but omit the results.
 \begin{table}[H]
 \caption{LFA predicted two-grid convergence factor, $\rho_h$,  as a  function of   $\nu$  and $h$  with $\beta=10^{-2}$ for exact stiffness-based BSR.}
\centering
\begin{tabular}{ l  c c c }
\hline
$h$                   & $\nu =2$    & $\nu=3$          &$\nu=4$      \\
\hline
$1/64$               &0.111         &0.037               &0.012         \\
$1/128$             &0.111         &0.037               &0.012         \\
$1/256$             &0.111         &0.037               &0.012        \\
\hline
\end{tabular}\label{tab:rho-nu}
\end{table}

\subsection{Multigrid performance}
In this subsection, we report some numerical results to confirm our theoretical results. We consider the distributed control problem on the unit square domain with homogeneous Dirichlet boundary conditions and $\hat{u}=0$ such that the discrete solution is  $z_h=0$.   The coarsest grid is a $4\times 4$ mesh.   For the stiffness matrix $A_{fd}$, we consider rediscretization operator on the coarse meshgrid.   We measure the multigrid convergence factor by the following form \cite{MR1807961}
\begin{equation*}
\hat{\rho}_{h} = \sqrt[n]{\frac{||d_h^{(n)}||_2}{||d_h^{(0)}||_2}},
\end{equation*}
where $d_h^{(n)}=b_h-\mathcal{L}_h z_h^{(n)}$, and $z_h^{(n)}$ is the approximation to the solution of \eqref{eq:saddle-structure}  at the
$n^{\rm th}$ multigrid iteration. We display the convergence factors obtained after 100 multigrid
cycles for different values of  $\nu$.   The initial guess
$z_h^{(0)}$ is chosen randomly.   We use the MATLAB code of Rees \cite{Rees2010github} to
generate the linear systems.

Table \ref{tab:BSR-dirichlet-beta4} shows the measured $W$-cycles multigrid convergence factors for exact  stiffness-based BSR  with algorithmic parameters,  $ \alpha_B=1, \omega_B=\frac{3}{4}$, for different values of $\beta$ and grid size $h$.  We see that our numerical results are independent of grid size $h$, and   $\hat{\rho}_h$ for $\nu>1$  matches well with the LFA two-grid convergence factor  $\rho_h$. Although for $\nu>1$ the measured convergence factors are  slightly better than the predictions, it is not surprising since LFA two-grid prediction is the same as the actual convergence for problems with periodic boundary conditions, and here we consider Dirichlet boundary conditions.

 \begin{table}[H]
 \caption{ $W$-cycles measured multigrid convergence factors, $\hat{\rho}_h$, for exact stiffness-based BSR with algorithmic parameters, $ \alpha_B=1, \omega_B=\frac{3}{4}$, compared with LFA two-grid predictions, $\rho_h$, for exact BSR with $ \alpha_B=1, \omega_B=\frac{3}{4}$.}
\centering
\begin{tabular}{ l  c c c c }
\hline
$\nu$  & 1    & 2   &3    & 4  \\
\hline
$\rho_h$            &0.333      &0.111         &0.037               &0.012       \\
\hline
\hline
    \multicolumn{5} {c}  {$\beta=10^{-2}$ } \\
\hline
$\hat{\rho}_{h=1/64}$              &0.279       &0.093         &0.031       &0.010         \\
$\hat{\rho}_{h=1/128}$            &0.275        &0.091         &0.030       &0.010           \\
$\hat{\rho}_{ h=1/256}$           &0.272       &0.090         & 0.030      & 0.010         \\
\hline
\hline
    \multicolumn{5} {c}  { $\beta=10^{-4}$ } \\
\hline
$\hat{\rho}_{h=1/64}$              &0.268       &0.090         &0.030       &0.010          \\
$\hat{\rho}_{h=1/128}$            &0.264       &0.089          &0.029        &0.010          \\
 $\hat{\rho}_{h=1/256}$            &0.261       &0.088         &0.029       &0.010          \\
\hline
\hline
  \multicolumn{5} {c}  { $\beta=10^{-6}$ } \\
\hline
$\hat{\rho}_{h=1/64}$              &0.289        & 0.094        &0.031       & 0.012        \\
$\hat{\rho}_{h=1/128}$            &0.276        &0.093         &0.030       &0.010          \\
$\hat{\rho}_{ h=1/256}$           &0.272       &0.092         &0.030       & 0.010         \\
\hline
\hline
  \multicolumn{5} {c}  { $\beta=10^{-8}$ } \\
\hline
$\hat{\rho}_{h=1/64}$              &0.302       & 0.098        &0.032       &0.011          \\
$\hat{\rho}_{h=1/128}$            &0.299       &0.098        &0.032       & 0.012         \\
$\hat{\rho}_{ h=1/256}$           &0.290       &0.096         &0.030       &0.012          \\
\hline
\end{tabular}\label{tab:BSR-dirichlet-beta4}
\end{table}

As mentioned before, using exact BSR, we need to solve the Schur complement system \eqref{eq:solution-of-precondtion-s1} exactly. In practice, an inexact solve is preferred.
First, we tried two or three sweeps of weighted Jacobi iteration and Gauss-Seidel iteration to approximate the solution of  \eqref{eq:solution-of-precondtion-s1}.  However, experiments showed that these approaches maintain the performance observed for exact BSR only for $\nu=1$, and there is a degradation for $\nu>1$. Moreover, these approaches are very sensitive to algorithmic parameters.   Instead, we consider solving the Schur
complement system by applying a multigrid $V(2,2 )$-cycle using weighted Jacobi relaxation with weight $\omega_J$. In this test, we  use  $\frac{\omega_B}{\alpha_B}=\frac{3}{4}$ because $1-\frac{\omega_B}{\alpha_B}\frac{a}{\hat{a}}$ is an eigenvalue of the relaxation error-propagation operator for both exact and inexact version. We  found that  there are many choices of $\alpha_B, \omega_B=\frac{3}{4}\alpha_B$ and $ \omega_J$  resulting in good multigrid results.  In our report, we consider $\alpha_B=1.5, \omega_B=\frac{3}{4}\alpha_B$ and $ \omega_J=0.8$.    We observe that using only  2 or 3 $V(2,2)$-cycles on the approximate Schur complement achieves convergence factors essentially matching those in Table \ref{tab:BSR-dirichlet-beta4} for $\nu=1,2$.  We display the inexact BSR results  in Table \ref{tab:IBSR-dirichlet-beta0}, where 3 $V(2,2)$-cycles are used. We see there is a slight degradation  of the measured convergence factors for inexact BSR for $\nu=3,4$ compared with the results in Table \ref{tab:BSR-dirichlet-beta4}. However, for $\nu=1,2$, the measured convergence factors $\hat{\rho}$ agree very well with the results in Table \ref{tab:BSR-dirichlet-beta4}, and are $h$-independent and robust to $\beta$. In practice, we recommand $W(1,0)$-cycle and $W(1,1)$-cycle for inexact BSR.

 \begin{table}[H]
 \caption{$W$-cycles measured convergence factors, $\hat{\rho}_h$,  for inexact stiffness-based BSR with inner 3 $V(2,2)$-cycles and algorithmic parameters, $\alpha_B=1.5, \omega_B=\frac{3}{4}\alpha_B, \omega_J=0.8$.}
\centering
\begin{tabular}{ l  c c c c }
\hline
$\nu$  & 1    & 2   &3    & 4  \\
\hline
\hline
    \multicolumn{5} {c}  {$\beta=10^{-2}$ } \\
\hline

$\hat{\rho}_{h=1/64}$              &0.280        &0.092         &0.046          &  0.044    \\
$\hat{\rho}_{h=1/128}$            &0.275         &0.092         &0.046         & 0.044          \\
$\hat{\rho}_{ h=1/256}$          &0.272         &0.090          &0.047        &0.045       \\
\hline
\hline
    \multicolumn{5} {c}  { $\beta=10^{-4}$ } \\
\hline
$\hat{\rho}_{h=1/64}$              &0.271         &0.089         &0.045        &  0.043          \\
$\hat{\rho}_{h=1/128}$            &0.267         &0.088         &0.046         & 0.044          \\
 $\hat{\rho}_{h=1/256}$           &0.264         &0.086         &0.047        &0.045          \\
\hline
\hline
  \multicolumn{5} {c}  { $\beta=10^{-6}$ } \\
\hline
$\hat{\rho}_{h=1/64}$              &0.283         &0.092          &0.039        &  0.025          \\
$\hat{\rho}_{h=1/128}$            &0.280       &0.091            &0.044        &  0.041           \\
$\hat{\rho}_{ h=1/256}$           &0.276        &0.090           &0.047        &0.045          \\
\hline
\hline
  \multicolumn{5} {c}  { $\beta=10^{-8}$ } \\
\hline
$\hat{\rho}_{h=1/64}$              & 0.296        &0.096         &0.032        &  0.010          \\
$\hat{\rho}_{h=1/128}$           &0.293          &0.095         &0.034         &0.025           \\
$\hat{\rho}_{ h=1/256}$           &0.289         &0.094         &0.043        &0.030          \\
\hline
\end{tabular}\label{tab:IBSR-dirichlet-beta0}
\end{table}

 Next, we consider the model problem \eqref{eq:Control-problem}  subject to the state equation \cite{rees2010optimal}
\begin{equation}\label{eq:sub-state-equ-nonzero}
\left\{
  \begin{aligned}
  -\triangle u &=&f,\,\,\,\, {\rm in}\,\,\Omega,\\
  u&=&\hat{u}|_{\partial \Omega},\,\,\, {\rm at}\,\,\, \partial \Omega,\\
   \end{aligned}
               \right.  \\
\end{equation}
where $\Omega=[0,1]^2$,  and
 \begin{equation*}
\hat{u} = \left\{
  \begin{aligned}
  & (2x-1)^2(2y-1)^2 & \quad   \text{if}\quad (x,y)\in[0,1/2]^2,\\
  & 0 &  \quad  \text{otherwise}.\\
   \end{aligned}
               \right.  \\
\end{equation*}
For this problem, we consider $W$-cycles multigrid methods  with inner 3 $V(2,2)$-cycles  on the Schur complement system with $\alpha_B=1.5, \omega_B=\frac{3}{4}\alpha_B, \omega_J=0.8$.  In Tables \ref{tab:Iteration-number-h-beta-nu1} and   \ref{tab:Iteration-number-h-beta-nu2}, we report the smallest number of iterations $n$   such that $\frac{||d_h^{(n)}||_2}{||d_h^{(0)}||_2}<10^{-10}$ for this distributed control model problem for different values of $\beta$, with $\nu=1$ and $\nu=2$, respectively. Again, from Tables \ref{tab:Iteration-number-h-beta-nu1} and   \ref{tab:Iteration-number-h-beta-nu2}, we see the robustness of our multigrid methods  with respect to the regularization parameters, $\beta$, and grid size $h$.

 \begin{table}[H]
 \caption{Iteration counts for the control model problem with  state equation \eqref{eq:sub-state-equ-nonzero} using  $\nu=1$ and $W$-cycles multigrid for inexact stiffness-based BSR with inner 3 $V(2,2)$-cycles and    algorithmic parameters, $\alpha_B=1.5, \omega_B=\frac{3}{4}\alpha_B, \omega_J=0.8$. }
\centering
\begin{tabular}{ l    c c c c}
\hline
$h$      & $\beta=10^{-2}$  &$\beta=10^{-4}$    &$\beta=10^{-6}$   & $\beta=10^{-8}$  \\
\hline
 $1/64$                     &17          & 17     &   17       &  19    \\
 $1/128$                   &  17             &17        &  17        &   18      \\
 $1/256$                   &18        &18      &18       &18   \\
  $1/512$                 &19       &19      & 19   &19\\
\hline
\end{tabular}\label{tab:Iteration-number-h-beta-nu1}
\end{table}

 \begin{table}[H]
 \caption{Iteration counts for the control model problem with  state equation \eqref{eq:sub-state-equ-nonzero}  using $\nu=2$ and $W$-cycles multigrid for inexact stiffness-based BSR with inner 3 $V(2,2)$-cycles and   algorithmic parameters, $\alpha_B=1.5, \omega_B=\frac{3}{4}\alpha_B, \omega_J=0.8$. }
\centering
\begin{tabular}{ l  c c c c c}
\hline
$h$      & $\beta=10^{-2}$  &$\beta=10^{-4}$    &$\beta=10^{-6}$   & $\beta=10^{-8}$  \\
\hline
 $1/64$                     & 12      & 12          &  11        &  10     \\
 $1/128$                   & 12        &  12      & 12          &  11      \\
 $1/256$                   &13         &13         & 12          &12    \\
  $1/512$                 &13         &13      & 13    &13 \\
\hline
\end{tabular}\label{tab:Iteration-number-h-beta-nu2}
\end{table}

\section{Conclusion} \label{sec:concl}
In this work, we present a novel multigrid method for the elliptic distributed control problems,  where stiffness-based  Braess-Sarazin multigrid relaxation scheme is proposed to the discrete saddle-point system of the control problems.  In this relaxation scheme, the stiffness matrix obtained from the five-point finite difference method for the Laplacian is used to approximate the mass matrix in the saddle-point system discretized by a finite element method. We hire LFA to minimize LFA smoothing factor for Braess-Sarazin  relaxation, leading to  a highly satisfactory convergence factor.   For practical use,  we develop  an inexact version of Braess-Sarazin  relaxation,  where a 3 $V(2,2)$-cycles multigrid with weighted Jacobi relaxation is applied to the Schur complement system. This inexact version preserves the optimal smoothing factor observed for exact Braess-Sarazin  relaxation.  Numerical results show the robustness of our multigrid methods with respect to the regularization parameters and grid size.  In future, we will extend the relaxation scheme considered here to  more complex optimal problems, such as the Stokes control problem.

\bibliographystyle{siam}