\section{Introduction}
Our aim in this paper is to design a simple consistent penalty method for
contact problems that avoids the solution
of variational inequalities. We eliminate the need for Lagrange multipliers to
enforce the contact conditions by using its definition in a discrete setting,
following an idea of Chouly and Hild \cite{CH13b} used for elastic contact.

\subsection{The model problem}

We consider the obstacle problem of finding the displacement $u$ of a membrane
constrained to stay above an obstacle given by $\psi=\psi(x,y)$ (with $\psi \leq 0$ at $\partial\Omega$):
\begin{equation}\label{obstacle}
\begin{array}{rcl}
-\Delta u-f &\geq& 0 \mbox{ in } \Omega\subset \mathbb{R}^2 \\
u &\geq&  \psi \mbox{ in } \Omega\\
(u-\psi) (f+\Delta u) &= & 0  \mbox{ in }  \Omega\\
u &= & 0  \mbox{ on } \partial \Omega
\end{array}
\end{equation}
where $\Omega$ is a convex polygon.
It is well known that this problem admits a unique solutions $u \in
H^1(\Omega)$. This follows from the theory of Stampacchia applied to
the corresponding variational inequality (see for instance
\cite{HHN96}).

\subsection{The finite element method}

There exists a large body of literature treating finite
element methods for unilateral problems in general and
obstacle problems in particular, e.g., \cite{KO88,Sch84,GleT89,Joh92,HHN96,ChNo00,Ve01,Zh01,BrCaHo07,WeWo10}.
Discretization of \eqref{obstacle} is usually performed directly starting
from the variational inequality or using a penalty method. The first
approach however leads to some nontrivial
choices in the construction of the discretization spaces in order to
satisfy the nonpenetration condition and associated inf-sup conditions
and until recently it has proved difficult to
obtain optimal error estimates \cite{HR12,DH16}. The latter approach, o
n the other hand leads to the usual
consistency and conditioning problems of penalty methods.

An alternative is to use the augmented Lagrangian method. We
introduce
the Lagrange multiplier $\lambda$ such that
\begin{equation}\label{mult1}
\begin{array}{rcl}
-\Delta u +\lambda &=& f \mbox{ in } \Omega \\
u &= & 0  \mbox{ on } \partial \Omega
\end{array}
\end{equation}
under the Kuhn-Tucker side conditions
\begin{equation}\label{mult2}
\begin{array}{rcl}
\psi -u &\leq&  0 \mbox{ in } \Omega\\
\lambda &\leq&  0 \mbox{ in } \Omega\\
(\psi-u) \,\lambda &= & 0  \mbox{ on }  \Omega .
\end{array}
\end{equation}
Using the standard trick of rewriting the Kuhn-Tucker conditions as
\begin{equation}
\lambda = -\frac{1}{\gamma}\left[ \psi - u   -\gamma\lambda\right]_+
\end{equation}
where $[x]_\pm = \pm \max(0,\pm x)$ and $\gamma\in \mathbb{R}^+$, cf., e.g., Chouly and Hild \cite{CH13b}, we can formulate the augmented Lagrangian problem of finding $(u,\lambda)$ that are stationary points to the functional
\begin{align}\nonumber
\frak{F}(u,\lambda) := {}& \frac12\int_{\Omega} \vert\nabla u\vert^2\, d\Omega + \int_{\Omega}\frac{1}{2\gamma}\left[  \psi-u-\gamma\lambda\right]_+^2\, d\Omega \\
 {}& -\frac12\int_{\Omega} \gamma \lambda^2 \, d{\Omega}
-\int_{\Omega} f u \, d{\Omega} ,
\end{align}
cf. Alart and Curnier \cite{AC91}, leading to seeking $(u,\lambda)\in H_0^1(\Omega)\times L_2(\Omega)$ such that
\begin{equation}
 \int_{\Omega} \nabla u\cdot\nabla v\, d\Omega - \int_{\Omega}\frac{1}{\gamma}\left[ \psi - u -\gamma\lambda\right]_+ v\, d\Omega =\int_{\Omega} f v \, d{\Omega}\quad\forall v\in H_0^1(\Omega)
\end{equation}
and
\begin{equation}
\int_{\Omega}\frac{1}{\gamma}\left[ \psi - u -\gamma\lambda\right]_+ \mu\, d\Omega +\int_{\Omega} \lambda \mu  \, d{\Omega}=0\quad\forall \mu\in L_2(\Omega) .
\end{equation}

For our discrete method, we assume that $\{ \mathcal{T}\}_h$ is a
family of conforming shape regular meshes on $\Omega$, consisting of
triangles $\mathcal{T} = \{T\}$
and
define $V_h$ as the space of $H^1$--conforming piecewise polynomial
functions
on $\mathcal{T}$, satisfying the homogeneous boundary condition of
$\Gamma_D$.
\[
V_h := \{v_h \in H^1_0(\Omega): v_h\vert_T \in \mathbb{P}_k(T), \,
\forall T \in \mathcal{T} \},\quad \mbox{ for } k \ge 2.
\]
We then formally replace $\lambda$ element--wise by
$\Delta u_h + f$ to obtain a discrete minimization problem:
seek $u_h\in V_h$ such that
\begin{equation}
u_h = \arg\min_{v\in V_h} \frak{F}_h(v)
\end{equation}
where
\begin{align}\nonumber
\frak{F}_h(v) := {}& \frac12\int_{\Omega} \vert\nabla v\vert^2\, d\Omega + \sum_{T\in\mathcal{T}}\int_{T}\frac{1}{2\gamma}\left[ \psi - v -\gamma(\Delta v + f)\right]_+^2\, d\Omega  \\
{}& -\frac12\sum_{T\in\mathcal{T}}\int_{T} \gamma (\Delta v + f)^2 \, d{\Omega}
-\int_{\Omega} f v \, d{\Omega} .\label{perturbed}
\end{align}
The
Euler--Lagrange equations corresponding to (\ref{perturbed}) take the form: Find $u_h \in V_h$ such that
\begin{equation}\label{FEM}
a(u_h,v_h) + b(u_h,\psi,f;v_h) =(f,v_h)_\Omega \quad \forall v_h
\in V_h
\end{equation}
where $(\cdot,\cdot)_\Omega$ denotes the standard $L^2$-inner product,
$a(u_h,v_h) := (\nabla u_h,\nabla v_h)_{\Omega}$ and
\begin{align}\nonumber
b(u_h,\psi,f;v_h):={}& \left<
 \gamma^{-1}[\psi -u_h- \gamma (\Delta u_h + f)]_+, -v_h - \gamma \Delta
 v_h\right>_{h} \\
 &{} -   \left<\gamma (\Delta u_h + f),  \Delta
 v_h\right>_{h} \label{stab_form0}
\end{align}
where
\[
\left<x_h,y_h\right>_{h} := \sum_{T\in\mathcal{T}} \int_T x_h y_h~ \mbox{d} x
\]
and, for use below,
\[
\| x_h\|_h := \left<x_h,x_h\right>_{h}^{1/2}.
\]
To simplify the notation below we introduce $P_\gamma(u_h) =
\gamma \Delta u_h+u_h $ and
\[
b(u_h,\psi,f;v_h):=\left<
 \gamma^{-1}[\psi - \gamma f-P_\gamma(u_h)]_+, -P_\gamma(v_h) \right>_{h} -   \left<\gamma (\Delta u_h + f),  \Delta
 v_h\right>_{h}.
\]
We will also omit $\psi$ and $f$ from the
argument of $b$ below, and use the notation $\Psi:=\psi -
\gamma f$ so that
\[
b(u_h;v_h):=\left<
 \gamma^{-1}[\Psi-P_\gamma(u_h)]_+, -P_\gamma(v_h) \right>_{h} -   \left<\gamma (\Delta u_h + f),  \Delta
 v_h\right>_{h}.
\]
Note that the form $b$ can be interpreted as a nonlinear consistent
least squares penalty term for the imposition of the contact
condition. A similar method was proposed in Stenberg et al. \cite{GuStVi16}
in the framework of variational inequalities.

We will below alternatively use the compact notation
\[
A_h(u_h,v_h):=a(u_h,v_h) + b(u_h;v_h)
\]
and the associated formulation, find $u_h\in V_h$ such that
\begin{equation}\label{CompFEM}
A_h(u_h;v_h) = (f,v_h)_\Omega,\mbox{ for all }
v_h\in V_h.
\end{equation}

\subsection{Summary of main results and outline} In Section 2 we
recall some technical
results, in Section 3 we derive an existence result for the discrete solution
using Brouwer's fixed point theorem and we prove uniqueness of the solution
using monotinicity of the the nonlinearity, in Section 4 we prove an a
priori error estimate and using a saturation assumption we also prove an
a posteriori error estimate, finally in Section 5 we present numerical
results confirming our theoretical
results and illustrating the performance of an adaptive algorithm
based on our a posteriori error estimate.

\section{Technical results}
Below we will use the notation $a \lesssim b$ for $a \leq C b$ where
$C$ is a constant independent of $h$, but not of the local mesh
geometry.

First we recall the following inverse inequality,
\begin{equation}\label{eq:inv_ineq}
\|\nabla v_h\|_{T} \leq C_i h_T^{-1} \|v_h\|_T , \qquad T \in \mathcal{T}
\end{equation}
see Thom\'ee \cite{Thomee2006}.

We will use the Scott-Zhang interpolant preserving boundary
conditions, denoted $i_h:H^1_0(\Omega) \mapsto V_h$. This operator is
$H^1$-stable, $\|i_h u\|_{H^1(\Omega)} \lesssim \|u\|_{H^1(\Omega)}$
and the following
interpolation error estimate is known to hold \cite{SZ90},
\begin{equation}\label{interpolation}
\|u-i_hu\|_{\Omega}+ h \|u-i_hu\|_{H^1(\Omega)}
+  h^2 \|
\Delta (u -i_hu)\|_{h} \lesssim  h^{k+1} |u|_{H^{k+1}(\Omega)}.
\end{equation}

The essential properties of the nonlinearity are collected in the
following lemmas.
\begin{lem}\label{lem:monotone}
Let $a,b \in \mathbb{R}$ then there holds
\[
([a]_+-[b]_+)^2 \leq ([a]_+ - [b]_+) (a - b),
\]
\[
|[a]_+-[b]_+| \leq |a-b|.
\]
\end{lem}
\begin{pf}
Developing the left hand side of the expression we have
\[
[a]_+^2 +[b]_+^2 - 2 [a]_+[b]_+ \leq  [a]_+ a + [b]_+ b -  a [b]_+
-  [a]_+b =  ([a]_+ - [b]_+) (a - b).
\]
For the proof of the second claim, this is trivially true in case both $a$
and $b$ are positive or negative. If $a$ is negative and $b$ positive then
\[
|[a]_+-[b]_+| = |b| \leq |b-a|
\]
and similarly if $b$ is negative and $a$ positive
\[
|[a]_+-[b]_+| = |a| \leq |b-a|.
\]
\end{pf}
\begin{lem}(Continuity of $b(\cdot;\cdot)$)\label{bcont}
For all $u_1,u_2,v \in V_h$, the form \eqref{stab_form0} satisfies
\begin{equation}
|b(u_1;v) - b(u_2;v)| \lesssim
\gamma^{-1} (\| (u_1-u_2)\|_{\Omega}+ \gamma h^{-1} \|\nabla (u_1-u_2)\|_\Omega) (\|v\|_{\Omega}+ \gamma h^{-1}
\|\nabla v_h\|_\Omega).
\end{equation}
\end{lem}
\begin{pf}
\begin{align*}
b(u_1;v_h) - b(u_2;v_h) = {}& \gamma^{-1}\left<
 [ \Psi-P_\gamma(u_1)]_+- [\Psi-P_\gamma(u_2)]_+,- P_\gamma(v_h) \right>_{h}  \\
{}& -   \left<\gamma \Delta (u_1 - u_2),  \Delta
 v_h\right>_{h}.
\end{align*}
Using the second inequality of Lemma \ref{lem:monotone} we see that
the nonlinearity satisfies
\begin{multline}\label{1st_ineq}
\gamma^{-1}|\left<
 [\Psi-P_\gamma(u_1) ]_+- [\Psi-P_\gamma(u_2)]_+, -P_\gamma(v_h)\right>_{h}| \\
\leq \gamma^{-1} \| P_\gamma(u_1) + \Psi- P_\gamma(u_2) -\Psi\|_{C,f} \|P_\gamma(v_h)\|_{h}.
\end{multline}
By the inverse inequality \eqref{eq:inv_ineq} we have
\begin{equation}\label{2nd_ineq}
\left<\gamma \Delta (u_1 - u_2),  \Delta
 v_h\right>_{h} \leq C_i^2 \gamma h^{-2} \|\nabla (u_1 - u_2)\|_\Omega
\|\nabla v_h\|_\Omega
\end{equation}
and
\begin{equation}\label{3rd_ineq}
\|P_\gamma(v_h)\|_{h} \leq \|v_h\|_\Omega + C_i \gamma h^{-1} \|\nabla v_h\|_\Omega.
\end{equation}
Collecting \eqref{1st_ineq},\eqref{2nd_ineq} and \eqref{3rd_ineq} we have
\begin{align*}
|b(u_1;v_h) - b(u_2;v_h) | \leq {}& (\|u_1 - u_2\|_\Omega + C_i \gamma
h^{-1} \|\nabla (u_1 - u_2)\|_\Omega) (\|v_h\|_\Omega + C_i \gamma
h^{-1} \|\nabla v_h\|_\Omega)\\
{}& + C_i^2 \gamma h^{-2} \|\nabla (u_1 - u_2)\|_\Omega
\|\nabla v_h\|_\Omega
\end{align*}
and the claim follows.
\end{pf}
\section{Existence of unique discrete solution}
In the previous works on Nitsche's method existence and uniqueness has
been proven by using the monotonicity and hemi-continuity of the
operator. Here we propose a different approach where we use the
Brouwer's fixed point theorem to establish existence and the
monotonicity of the nonlinearity for uniqueness. We start by showing
some positivity results and a priori bounds. Since we are interested
in existence and uniqueness for a fixed mesh parameter $h$, we do not
require that the bounds in this section are uniform in $h$.
\begin{lem}\label{lem:monotonicity}
Let $u_1,u_2 \in V_h$ and assume that
\begin{equation}\label{pos_cond}
\gamma < C_i^{-2} h^2/2
\end{equation}
 then there holds
\begin{multline*}
\frac{\alpha}{2} \|u_1-u_2\|_{H^1(\Omega)}^2 +\gamma^{-1} \|[\Psi-P_\gamma(u_1) ]_+- [\Psi-P_\gamma(u_2)]_+\|^2_{h}\\
\leq A_h(u_1;u_1-u_2)-A_h(u_2;u_1-u_2)
\end{multline*}
and
\[
\frac{\alpha}{4} \|u_1\|_{H^1(\Omega)}^2
\leq A_h(u_1;u_1) + C \alpha^{-1} \gamma^{-2} \|[\Psi]_+\|^2_\Omega.
\]
\end{lem}
\begin{pf}
First observe consider the form $b(\cdot;\cdot)$,
\begin{multline*}
b(u_1;v_h) - b(u_2;v_h) \\=\gamma^{-1}\left<
 [ \Psi-P_\gamma(u_1)]_+- [\Psi-P_\gamma(u_2)]_+, -v_h-\gamma \Delta
 v_h + \Psi - \Psi\right>_{h}  \\
-   \left<\gamma \Delta (u_1 - u_2),  \Delta
 v_h\right>_{h}.
\end{multline*}
Using the monotonicity of Lemma \ref{lem:monotone}  we may write
\begin{align*}
b(u_1;u_1-u_2) - b(u_2;u_1-u_2)
\ge {}& \gamma^{-1}\|[ \Psi-P_\gamma(u_1) ]_+- [\Psi - P_\gamma(u_2) ]_+\|^2_{h}\\ &{} -   \gamma \|\Delta (u_1 - u_2)\|^2_{h}.
\end{align*}
Observe that using an inverse inequality \eqref{eq:inv_ineq} we have $$\gamma \|\Delta (u_1 - u_2)\|^2_{h} \leq C_i \gamma
h^{-2} \|\nabla(u_1 - u_2) \|^2_\Omega$$.
We may then write
\begin{multline*}
(1 -  C^2_i
h^{-2}\gamma) \|\nabla (u_1-u_2)\|_{\Omega}^2 + \gamma^{-1}\|[\Psi-P_\gamma(u_1)
]_+- [\Psi-P_\gamma(u_2) ]_+\|^2_{h}\\
\leq A_h(u_1;u_1-u_2)-A_h(u_2;u_1-u_2)
\end{multline*}
It follows that choosing
$\gamma < C_i^{-2} h^2 2$ and applying the Poincar\'e inequality
\begin{equation}\label{eq:poincare}
\alpha^{\frac12} \|u\|_{H^1(\Omega)} \leq \|\nabla u\|_\Omega, \quad
\forall u \in H^1_0(\Omega)
\end{equation}
 there holds
\begin{multline*}
\frac{\alpha}{2} \|u_1-u_2\|_{H^1(\Omega)}^2 + \gamma^{-1}\|[\Psi-P_\gamma(u_1) ]_+- [\Psi-P_\gamma(u_2) ]_+\|^2_{h}\\
\leq A_h(u_1;u_1-u_2)- A_h(u_2;u_1-u_2).
\end{multline*}
The second inequality follows by taking $u_2 = 0$ above and noting that
then
\begin{align*}
&\frac{\alpha}{2} \|u_1\|_{H^1(\Omega)}^2 +\gamma^{-1} \|[\Psi-P_\gamma(u_1)]_+- [\Psi]_+\|^2_{h}
\\
&\qquad \leq A_h(u_1;u_1)- \gamma^{-1} \left< [\Psi]_+,- P_\gamma(
 u_1) \right>_{h} \\
&\qquad \leq A_h(u_1;u_1) + (1+C_i \gamma h^{-1})\alpha^{-\frac12} \gamma^{-1}
\|[\Psi]_+\|_\Omega  \alpha^{\frac12} \|u_1\|_{H^1(\Omega)}
\end{align*}
where we used \eqref{3rd_ineq} in the last step. Considering the
condition on $\gamma$ and using an arithmetic-geometric inequality we
may conclude.
\end{pf}
\begin{prop}
Assume that  $\gamma$ saisfy \eqref{pos_cond}.
Then formulation \eqref{CompFEM} using the contact operator
\eqref{stab_form0}, admits a unique solution.
\end{prop}
\begin{pf}
The uniqueness is an immediate consequence of Lemma \ref{lem:monotonicity}. If $u_1$ and $u_2$ both are solution to \eqref{CompFEM},
then
\[
A_h(u_1;u_1-u_2)-A_h(u_2;u_1-u_2) = (f,u_1-u_2)_\Omega
-(f,u_1-u_2)_\Omega = 0
\]
and we conclude that $\|u_1-u_2\|_{H^1(\Omega)}=0$ and hence
$u_1\equiv u_2$.
Let $N_V$ denote the number of
degrees of freedom of $V_h$.

Consider the mapping $G:\mathbb{R}^{N_V} \mapsto \mathbb{R}^{N_V}$ defined by
\[
(G(U),V)_{\mathbb{R}^{N_V}}  := A_h(u_h,v_h)-
(f,v_h)_\Omega,
\]
where $U = \{u_i\}_{i=1}^{N_V}$, $V=\{v_i\}_{i=1}^{N_V}$,
where $\{u_i\}$ and $\{v_i\}$ denotes the vectors of unknown
associated to the basis functions of $V_h$.

By the second claim of Lemma \ref{lem:monotonicity},
there holds
\[
\frac{\alpha}{4} \|u_h\|^2_{H^1(\Omega)} - C\alpha^{-1} \gamma^{-2}  \|[\Psi]_+\|_\Omega^2 -  (f,u_h)_\Omega \leq
A_h(u_h,u_h) - (f,u_h)_\Omega = (G(U),U)_{\mathbb{R}^{N_V}}.
\]
Since
\[
\frac{\alpha}{4} \|u_h\|^2_{H^1(\Omega)}  - (f,u_h)_\Omega \ge
\frac{\alpha}{8} \| u_h\|^2_{H^1(\Omega)} -  C\frac{1}{\alpha} \|f\|_\Omega^2
\]
we have that for any fixed $h$ the following positivity holds for $U$ sufficiently large
\[
0 < \frac{\alpha}{8} \|u_h\|^2_{H^1(\Omega)}  - \frac{C}{\alpha}
(\gamma^{-2}\|[\Psi]_+\|_\Omega^2+\|f\|_\Omega^2) \leq (G(U),U)_{\mathbb{R}^{N_V}}.
\]
Assume that this positivity holds whenever $|U|\ge q \in
\mathbb{R}_+$. Denote by $B_q$ the (closed) ball in
$\mathbb{R}^{N_V}$ with radius $q$
and assume that there is no $U \in B_q$ such that $G(U)=0$. Define the
function
\[
\phi(U) = - q G(U)/|G(U)|.
\]
Then $\phi:B_{q} \mapsto B_{q}$, $\phi$ is
continuous by Lemma \ref{bcont} and the assumption that $|G(U)|>0$ for
all $U \in B_q$. Hence there
exists a fixed point $X \in B_q$ such that
\[
X = \phi(X).
\]
It follows that
\[
|X|^2 = -q (G(X),X)/|G(X)|,
\]
but since $|X|=q$, by assumption $(G(X), X)>0$, which leads to a
contradiction, since $|X|>0$.
\end{pf}
\section{Error estimates}

\begin{thm}(A priori error estimate)\label{thm:apriori}
Assume that $u \in H^1_0(\Omega)$ with $\Delta u \in L^2(\Omega)$ is
the solution of \eqref{obstacle} and that $u_h$ is the solution to \eqref{FEM}
with \eqref{stab_form0} and $0<\gamma =\gamma_0 h^{2}$,
where $\gamma_0 \in \mathbb{R}$, $\gamma_0 < C_i^{-2}/2$.
then there holds for all $v_h \in
V_h $
\begin{multline}\label{best_approximation}
\alpha\|u-u_h\|_{H^1(\Omega)}^2 +\gamma^{-1} \|  [\Psi-P_{\gamma}(u_h)]_+- [\Psi-P_{\gamma}(u)]_+\|^2_{h}  \\
\lesssim \frac{1}{ \alpha} \|u-v_h\|_{H^1(\Omega)}^2+
\|\gamma^{-\frac12}(u-v_h)\|_{\Omega}
+  \|\gamma^{\frac12}
\Delta (u -  v_h)\|_{h}^2.
\end{multline}
If in addition $u \in H^{k+1}(\Omega)$ then there holds
\begin{equation}\label{optimal_convergence}
\alpha\|u-u_h\|_{H^1(\Omega)} +\gamma^{-1/2} \|
[\Psi-P_{\gamma}(u_h)]_+- [\Psi-P_{\gamma}(u)]_+\|_{h} \lesssim h^k |u|_{H^{k+1}(\Omega)}.
\end{equation}
\end{thm}
\begin{pf}
Using the definition of $a(\cdot,\cdot)$ we may write
\begin{align*}
\|\nabla (u-u_h)\|_{\Omega}^2
&\leq a(u-u_h,u-u_h)
\\
&= a(u-u_h,u - v_h) +a(u-u_h,v_h-u_h)
\\
&\leq \frac{\alpha}{4} \|u-u_h\|_{H^1(\Omega)}^2
+ \frac{1}{\alpha} \|u-v_h\|_{H^1(\Omega)}^2 +a(u-u_h,v_h-u_h).
\end{align*}
Observe that
\begin{align}\nonumber
a(u,v_h-u_h) = {}& \left<-\Delta u - f + f, v_h-u_h\right>_{V',V}
\\
=  {}& \left<\gamma^{-1} [\Psi-P_{\gamma}(u)]_+,
 (v_h-u_h) \right>_{V',V} + (f, v_h-u_h)_{\Omega}.
\end{align}
If $ [\Psi-P_{\gamma}(u)]_+ \in L^2(\Omega)$ we may also write
\[
\left<\Delta u + f, \gamma \Delta(v_h-u_h)\right>_h+\left<\gamma^{-1} [\Psi-P_{\gamma}(u)]_+,
 \gamma \Delta(v_h-u_h)\right>_h = 0.
\]
It follows that
\begin{align}\nonumber
a(u,v_h-u_h)
= {}& (f, v_h-u_h)_{\Omega} - \left<\gamma^{-1} [\Psi-P_{\gamma}(u)]_+,- P_{\gamma}(v_h-u_h)\right>_h \\ \nonumber &{} + \left<\Delta u + f, \gamma \Delta (v_h-u_h)\right>_h\\
= {}& (f,v_h-u_h)_{\Omega} - b(u;v_h-u_h).
\end{align}
As a consequence we have the following property reminiscent of Galerkin
orthogonality,
\begin{align}\nonumber
&a(u-u_h,v_h-u_h)
\\ \nonumber
&\qquad =  b(u_h;v_h-u_h)- b(u;v_h-u_h)
\\ \label{eq:GO_pert}
&\qquad =\left<\gamma^{-1}
  [\Psi-P_{\gamma}(u_h)]_+-\gamma^{-1} [\Psi-P_{\gamma}(u)]_+,
 -P_{\gamma}(v_h-u_h) \right>_{h}
 \\ \nonumber
&\qquad \qquad - \gamma \left<\Delta (u_h - u), \Delta (v_h - u_h) \right>_{h}
\end{align}
First observe that
\begin{align*}
&\gamma\left<\Delta (u_h - u), \Delta (v_h - u_h) \right>_{h}
\\
&\qquad
\leq
\|\gamma^{\frac12}(\Delta u_h - \Delta v_h )\|^2_{h} + \|\gamma^{\frac12}
(\Delta v_h - \Delta u )\|_{h}\|\gamma^{\frac12} (\Delta v_h - \Delta u_h )\|_{h}
\\
&\qquad
\leq
\frac12 C^2_i h^{-2} \gamma \|\nabla (u_h - v_h)\|^2_\Omega +  \|\gamma^{\frac12}
(\Delta v_h - \Delta u )\|_{h}^2
\\
&\qquad
\leq
C^2_i h^{-2} \gamma \|\nabla (u_h - u)\|^2_\Omega + C^2_i h^{-2}
\gamma \|\nabla (v_h - u)\|^2_\Omega + \|\gamma^{\frac12}
(\Delta v_h - \Delta u )\|_{h}^2.
\end{align*}
Considering the first term in the right hand side of equation
\eqref{eq:GO_pert} we may write
\begin{align*}
&\left<\gamma^{-1}
  [\Psi-P_{\gamma}(u_h)]_+-\gamma^{-1} [\Psi-P_{\gamma}(u)]_+,
  -P_{\gamma}(v_h-u_h) \right>_{h} \\
\\
&\qquad = \underbrace{\left<\gamma^{-1}
  [\Psi-P_{\gamma}(u_h)]_+-\gamma^{-1}
  [\Psi-P_{\gamma}(u)]_+,
  -P_{\gamma}(v_h-u) \right>_{h} }_{I}
\\
&\qquad \qquad + \underbrace{\left<\gamma^{-1}
  [\Psi-P_{\gamma}(u_h)]_+-\gamma^{-1}
  [\Psi-P_{\gamma}(u)]_+,
  -P_{\gamma}(u-u_h)  \right>_{h} }_{II}
\\
&\qquad = I+II.
\end{align*}
The term $I$ may be bounded using the Cauchy-Schwarz inequality
followed by the arithmetic geometric inequality
\[
I \leq \epsilon \gamma^{-1}\| [\Psi-P_{\gamma}(u_h)]_+- [\Psi-P_{\gamma}(u)]_+\|_{h}^2 + \frac{1}{4 \epsilon}
\|\gamma^{-\frac12} P_{\gamma}(v_h-u)\|_{h}^2.
\]
For the term $II$ we use the monotonicity property $([a]_+-[b]_+)(b-a)
\leq -([a]_+-[b]_+)^2$, with $a=\Psi-P_{\gamma}(u_h)$ and
$b=\Psi-P_{\gamma}(u)$ so that
\[
([a]_+-[b]_+)(b-a) = ([\Psi-P_{\gamma}(u_h)]_+-\gamma^{-1}
  [\Psi-P_{\gamma}(u)]_+)(
  \Psi -P_{\gamma}(u) - \Psi + P_\gamma(u_h))
\] to deduce that
\[
II \leq -\gamma^{-1} \|  [\Psi-P_{\gamma}(u_h)]_+- [\Psi-P_{\gamma}(u)]_+\|^2_{h}.
\]
Collecting the above bounds and using the Poincar\'e inequality
\eqref{eq:poincare} we find,
\begin{multline}
\alpha \left(\frac34 - C^2_i h^{-2} \gamma \right) \|u-u_h\|_{H^1(\Omega)}^2 \\+ \left(1-\epsilon \right) \gamma^{-1} \|  [\Psi-P_{\gamma}(u_h)]_+- [\Psi-P_{\gamma}(u)]_+\|^2_{h}
+ \left(1-\epsilon \right) \|\gamma^{\frac12}\Delta (u - u_h) \|^2_{h}  \\
\leq \frac{1}{\alpha} \|u-v_h\|_{H^1(\Omega)}^2+ \frac{1}{4 \epsilon}
\|\gamma^{-\frac12} P_{\gamma}(u-v_h)\|_{h}^2
+ \|\gamma^{\frac12} \Delta (u-v_h) \|^2_{h}
\end{multline}
Fixing $\epsilon=\frac12$, and fixing $\gamma$ sufficiently small so that $C_i^2h^{-2}
\gamma \leq \alpha/4$ then there holds
\begin{multline}
\alpha\|u-u_h\|_{H^1(\Omega)}^2 \gamma^{-1} \|  [\Psi-P_{\gamma}(u_h)]_+- [\Psi-P_{\gamma}(u)]_+\|^2_{h}  \\
\lesssim \frac{1}{ \alpha} \|u-v_h\|_{H^1(\Omega)}^2+
\|\gamma^{-\frac12}(u-v_h)\|_{h}^2
+  \|\gamma^{\frac12}
\Delta (u -v_h)\|_{h}^2.
\end{multline}
This concludes the proof of \eqref{best_approximation}. The error
estimate \eqref{optimal_convergence} then follows by choosing $v_h$ to
be the interpolant, $i_h u$, applying the approximation error
estimate \eqref{interpolation} on the form
\begin{align*}
&\|u-i_hu\|_{H^1(\Omega)}+
\|\gamma^{-\frac12}(u-i_hu)\|_{h}
+  \|\gamma^{\frac12}
\Delta (u -i_hu)\|_{h}
\\
&\qquad \lesssim ( h^{k}+ \gamma^{-1/2} h^{k+1}+
\gamma^{1/2} h^{k-1}) |u|_{H^{k+1}(\Omega)}
\end{align*}
and using the bound on $\gamma$.
\end{pf}
{\bf Assumption:} \emph{(Saturation) We assume that there exists a constant
$C_s$ such that}
\begin{equation}\label{saturation}
\|\Delta (u - u_h)\|_{h} \leq C_s h^{-1} \|\nabla (u - u_h)\|_\Omega.
\end{equation}
\begin{thm}(A posteriori error estimate)
Assume that $u \in H^1_0(\Omega)$ with $\Delta u \in L^2(\Omega)$ is
the solution of \eqref{obstacle} and $u_h$ the solution of \eqref{FEM}
satisfying \eqref{saturation} and with the parameter $\gamma$ satisfying $\gamma \leq \tfrac12 C_s^{-2} h^2$, then
\begin{equation}
\alpha\|u-u_h\|_{H^1(\Omega)}+\gamma^{-1} \|  [\Psi-P_{\gamma}(u_h)]_+- [\Psi-P_{\gamma}(u)]_+\|_{h} \lesssim E(h,\gamma,u_h,f),
\end{equation}
where
\[
E(h,\gamma,u_h,f):= h \|f+\Delta u_h +\gamma^{-1}
[\Psi-P_{\gamma}(u_h)]_+\|_{h} +  \|h^{\frac12} \jump{\partial_n u_h}\|_{\mathcal{F}}.
\]
\end{thm}
\begin{pf}
Let $e=u-u_h$ then, under the assumption \eqref{saturation} and using
\eqref{lem:monotonicity} we may write
\begin{align}\nonumber
&(1 - C_s^2 h^{-2}\gamma )\|\nabla e\|_\Omega^2 + \gamma^{-1} \|  [\Psi-P_{\gamma}(u)]_+- [\Psi-P_{\gamma}(u_h)]_+\|^2_{h}
\\ \label{eq:sat_bound1}
&\qquad \leq (\nabla (u-u_h),\nabla e)_\Omega
\\ \nonumber
&\qquad \qquad + \gamma^{-1} \left< [\Psi-P_{\gamma}(u)]_+- [\Psi-P_{\gamma}(u_h)]_+,- P_{\gamma}(e) \right>_{h}
\\ \nonumber
&\qquad \qquad -\gamma \left< \Delta (u - u_h), \Delta e\right>_{h}
\end{align}
Now, using similar arguments as in Theorem \ref{thm:apriori} we deduce
\[
 (\nabla u,\nabla e)_\Omega - \gamma^{-1}  \left<
   [\Psi-P_{\gamma}(u)]_+, P_{\gamma}(e) \right>_{h}-\gamma \left<
\Delta u, \Delta e \right>_{h} = \left<f,P_{\gamma}(e) \right>_{h}.
\]
Therefore the bound \eqref{eq:sat_bound1} may be written, under the assumption
$C_s^2 h^{-2} \gamma \leq 1/2$,
\begin{align}\nonumber
&\frac12 \|\nabla e\|_\Omega^2 + \gamma^{-1} \|  [\Psi-P_{\gamma}(u)]_+- [\Psi-P_{\gamma}(u_h)]_+\|^2_{h}
\\ \label{eq:sat_bound2}
&\qquad \leq
\left<f,P_{\gamma}(e) \right>_{h} - (\nabla u_h,\nabla e)_\Omega
\\ \nonumber
&\qquad\qquad  + \gamma^{-1}
\left< [\Psi-P_{\gamma}(u_h)]_+, P_{\gamma}(e) \right>_{h}+\gamma\left<
    \Delta u_h , \Delta e \right>_{h}
\\ \nonumber
&\qquad = \left<f+\Delta u_h +\gamma^{-1} [\Psi-P_{\gamma}(u_h)]_+,(I +\gamma \Delta) (e - i_h e) \right>_h
\\ \nonumber
&\qquad \qquad + \left<\jump{\partial_n u_h},e - i_h e
\right>_{\mathcal{F}}
\end{align}
Using that $\gamma \|\Delta i_h e\|_h \lesssim \gamma h^{-1} \|\nabla
i_h e\|_\Omega$, the $H^1$ stability of $i_h$ and the definition
of $\gamma$ we obtain
\begin{align*}
&\frac12 \|\nabla e\|_\Omega^2 + \gamma^{-1} \| [\Psi-P_{\gamma}(u)]_+- [\Psi-P_{\gamma}(u_h)]_+\|^2_{h}
\\
&\qquad
\leq
C (h +\gamma h^{-1}+\gamma^{\frac12})( \|f+\Delta u_h -\gamma^{-1} [\Psi-P_{\gamma}(u_h)]_+\|_{h}
\\
&\qquad \qquad \qquad \times (\|\nabla e\|_\Omega + \gamma^{\frac12} \|\Delta e\|_h)
\\
&\qquad \qquad+ C \|h^{\frac12} \jump{\partial_n u_h}\|_{\mathcal{F}} \|\nabla
e\|_\Omega
\\
&\qquad
\leq
C (h +\gamma h^{-1} +\gamma^{\frac12})( \|f+\Delta u_h
-\gamma^{-1} [\Psi-P_{\gamma}(u_h)\|_{h} + \|h^{\frac12}
\jump{\partial_n u_h}\|_{\mathcal{F}})
\\
&\qquad \qquad \qquad \times(1+
C_s \gamma^{\frac12} h^{-1}) \|\nabla
e\|_\Omega
\end{align*}

Once again applying the assumption on $\gamma$ we obtain the bound
\begin{multline*}
\|\nabla e\|_\Omega + \gamma^{-1} \| [\Psi-P_{\gamma}(u)]_+- [\Psi-P_{\gamma}(u_h)]_+\|_{h} \\
\leq  C h \|f+\Delta u_h -\gamma^{-1}
[\Psi-P_{\gamma}(u_h)]_+\|_{h} +  C\|h^{\frac12} \jump{\partial_n u_h}\|_{\mathcal{F}}.
\end{multline*}
We conclude the proof using the Poincar\'e inequality \eqref{eq:poincare}.
\end{pf}
\begin{rem}
This a posteriori error estimate has the disadvantage of the
saturation assumption and also that the parameter $\gamma$ depends on
the constant in the saturation assumption. However as we shall see
below it appears to give a very good representation of the $H^1$-error
and can be used to drive adaptive refinement.
\end{rem}

\section{Numerical examples}

\subsection{A smooth rotational symmetric exact solution}
This example, from \cite{NoSiVe03}, is posed on the square $\Omega=(-1,1)\times(-1,1)$ with $\psi=0$ and
\[
f=\left\{\begin{array}{c}
    -8 r_0^2 (1-(r^2-r_0^2)) \quad \text{if $r\leq r_0$},\\
  -8 (r^2+ (r^2-r_0^2)) \quad \text{if $r> r_0$},
\end{array}\right.
\]
where $r=\sqrt{x^2+y^2}$ and $r_0=1/4$, and
with Dirichlet boundary conditions taken from the corresponding exact solution
\[
u = [r^2-r_0^2]_+^2 .
\]

We choose $\gamma=\gamma_0 h^2$ with $\gamma_0 = 1/100$ and show the convergence in the $L_2$-- and $H^1$--norms in Figure \ref{fig:smooth} together with the
error indicator. We remark that the smoothness of the solution precludes mesh zoning for this example, but that the indicator has the same asymptotic behaviour as the $H^1$ error. An elevation of the computed solution on one of the meshes in a sequence is given in Fig. \ref{fig:elevetionsmooth}. We note the optimal convergence of $O(h^3)$ in $L_2$ and $O(h^2)$ in $H^1$. Here we use $h=1/\sqrt{\text{NNO}}$, wehre NNO denotes the number of nodes in a uniformly refined mesh.

\subsection{A non-smooth exact solution}

This example, from \cite{BrCaHo07}, is posed on the L-shaped domain $\Omega= (-2, 2)\times (-2, 2) \setminus [0, 2)\times (-2, 0]$ with $\psi=0$ and
\[
f(r,\varphi) = -r^{2/3}\sin{(2\varphi/3)}(\gamma'(r)/r+\gamma''(r))-\frac{4}{3}r^{-1/3}\gamma'(r)\sin(2\varphi/3)-\gamma_2(r)
\]
(note the sign error in \cite{BrCaHo07}), where, with $\hat{r}=2(r-1/4)$,
\[
\gamma_1(r)=\left\{\begin{array}{ll}1,& \hat{r} < 0\\
-6\hat{r}^5+15\hat{r}^4-10\hat{r}^3+1, & 0\leq\hat{r}<1\\
0, & \hat{r}\geq 1,\end{array}\right.
\]
\[
\gamma_2(r)=\left\{\begin{array}{ll}0,& r\leq 5/4 ,\\
1 & \text{elsewhere.}\end{array}\right.
\]
with Dirichlet boundary conditions taken from the corresponding exact solution
\[
u(r,\varphi)=r^{2/3}\gamma_1(r)\sin(2\varphi/3)
\]
which belongs to $H^{5/3-\varepsilon}(\Omega)$ for arbitrary $\varepsilon > 0$.

For this example we plot, in Fig. \ref{fig:nonsmootherror} the error on consecutive adaptively refined meshes, using the minimum meshsize as a measure of $h$.
We note the suboptimal convergence and that the indicator still approximately follows the $H^1$ error asymptotically. In Fig. \ref{fig:meshes}
we show the corresponding sequence of refined meshes, and in Fig. \ref{fig:elnonsmo} we show an elevation of the approximate solution on the final mesh in the sequance.

\bibliographystyle{plain}


\section{Introduction}
There have been theoretical approches to gravity theories, where a fundamental constraint is implemented. For example in the Two Measures Theories [1] one works, in addition to the regular measure of integration in the action $ \sqrt{-g} $, also with another measure which is also a denstiy and which is also a total derivative. In this case, one can use for constructing this measure 4 scalar fields $ \varphi_{a} $, where $ a=1,2,3,4 $. Then we can define the density $ \Phi=\varepsilon^{\alpha\beta\gamma\delta}\varepsilon_{abcd}\partial_{\alpha}\varphi_{a}\partial_{\beta}\varphi_{b}\partial_{\gamma}\varphi_{c}\partial_{\delta}\varphi_{d} $, and then we can write an action that uses both of these densities:

\begin{equation}
S=\int d^{4}x\Phi\mathcal{L}_{1}+\int d^{4}x\sqrt{-g}\mathcal{L}_{2}
\end{equation}
As a consequence of the variation with respect to the scalar fileds $ \varphi_{a} $, assuming that $ \mathcal{L}_{1} $
 and $ \mathcal{L}_{2} $
  are independed of the scalar fileds $\varphi_{a} $, we obtain that:

\begin{equation}
A_{a}^{\alpha}\partial_{\alpha}\mathcal{L}_{1}=0
\end{equation}

 where $ A_{a}^{\alpha}=\varepsilon^{\alpha\beta\gamma\delta}\varepsilon_{abcd}\partial_{\beta}\varphi_{b}\partial_{\gamma}\varphi_{c}\partial_{\delta}\varphi_{d} $. Since $ det[A_{a}^{\alpha}]\sim\Phi^{3} $ as one easily see, then that for $ \Phi\neq0 $ ,(2) implies that $ \mathcal{L}_{1}=M=Const $. This result can expresed as a covariant conservation of a stress energy momentum of the form $ T_{\left(\Phi\right)}^{\mu\nu}=\mathcal{L}_{1}g^{\mu\nu} $, and using the 2nd order formalism, where the covariant derivative of $ g_{\mu\nu} $ is zero, we obtain that $ \nabla_{\mu}T_{\left(\Phi\right)}^{\mu\nu}=0 $, implies ${\partial_{\alpha}\mathcal{L}_{1}=0}$.
   This suggest generalizing the idea of the Two Measures Theory, by imposing the covariant conservation of a more nontrivial kind of energy momenetum tensors, which we denote as $ T_{\left(\chi\right)}^{\mu\nu} $ [2].
   Therfore we consider an action of the form:
\begin{equation}
S=\mathcal{S}_{\left(\chi\right)}+\mathcal{S}_{\left(R\right)}=\int d^{4}x\sqrt{-g}\chi_{\mu;\nu}T_{\left(\chi\right)}^{\mu\nu}+\frac{1}{16\pi G}\int d^{4}x\sqrt{-g}R
\end{equation}

   where $ \chi_{\mu;\nu}=\partial_{\nu}\chi_{\mu}-\Gamma_{\mu\nu}^{\lambda}\chi_{\lambda} $. If we assume $ T_{\left(\chi\right)}^{\mu\nu} $ to be independent of$  \chi_{\mu} $
     and having $ \Gamma_{\mu\nu}^{\lambda} $ being defined as the Christoffel Connection Coefficients, then the variation with respect to $ \chi_{\mu} $ gives a covariant conservation $ \nabla_{\mu}T_{\left(\chi\right)}^{\mu\nu}=0 $. Notice the fact that the energy density is the canonically conjugated variable to $ \chi^{0} $, which is what we expect from a dynamical time (here represented by index of time).

\section{Solutions with homothetic vector filed}
Interesting particular cases of cosmological solutions are obtained when the vector filed $ \chi^{\mu} $  has the property of being a homothetic vector filed:

\begin{equation}
\mathcal{L}_{\chi} g_{\mu\nu} \equiv \chi_{\mu;\nu}+\chi_{\nu;\mu}=\frac{\chi}{2}g_{\mu\nu}
\end{equation}

when $ \chi=\chi_{;\lambda}^{\lambda} $ is a constant trace of $ \chi_{\mu;\nu} $ (for homothetic Killing vector). \\The effective stress energy momentum tensor, from variation of the metric using the action (3) (using the identities in appendix A) is:

\begin{equation}
 T_{eff}^{\mu\nu}=\frac{\chi T(\chi)}{8}g^{\mu\nu}+\frac{\chi}{4}g_{\alpha\beta}\frac{\partial}{\partial g_{\mu\nu}}T_{(\chi)}^{\alpha\beta}-\frac{1}{4}\chi T_{\left(\chi\right)}^{\mu\nu}-\frac{1}{2}\left(\chi^{\lambda}\nabla_{\lambda}\right)T_{(\chi)}^{\mu\nu}
\end{equation}

when $ T(\chi)=T_{\lambda}^{\lambda}(\chi) $  is the trace of $ T_{(\chi)}^{\mu\nu} $. We see that the stress-energy tensor $ T_{eff}^{\mu\nu} $  is related to original $ T_{(\chi)}^{\mu\nu} $ but it is not equal to it. Where $ T_{eff}^{\mu\nu} $  is the stress-energy momentum tensor that apear in the right side of Einstein equation  $T^{\mu\nu}_{eff}=R^{\mu\nu}-\frac{1}{2}g^{\mu\nu}R$. As we can see, there is a simple linear relation between $  \frac{1}{4}\chi T_{\left(\chi\right)}^{\mu\nu} $, and $ \frac{1}{2}\left(\chi^{\lambda}\nabla_{\lambda}\right)T_{(\chi)}^{\mu\nu} $  in some spaces the last expression could be proportional to $ T_{\left(\chi\right)}^{\mu\nu} $
 if such a tensor is an homogeneous function, by using Euler's homogeneity theorem, which it is satisfied for power law functions of time (in Cosmological solutions). In addition to conservation of $ T_{\left(\chi\right)}^{\mu\nu} $
  the original $ T_{eff}^{\mu\nu} $ also should be conserved. This two conditions are co-exist on the same metric of space-time.

In a context of Cosmology, because of the presence of the term $ \frac{1}{2}\left(\chi^{\lambda}\nabla_{\lambda}\right)T_{(\chi)}^{\mu\nu} $
  we are interested in using the Friedman-Lemaitre-Robertson-Walker (FLRW) metric, with spacial curvature $ k=0 $, wich is the only one that gives the property of proportionality between those two kinds of stress-energy tensors $ T_{\left(\chi\right)}^{\mu\nu} $ and $ T_{eff}^{\mu\nu} $. Otherwise, we get a contradiction between the two conservation of stress-energy momentum tensors. We find this mathematical condition in our universe, that from measurements we sure that there isn't spacial curvature, or if it's exists, it is very small. From our theory, this is a natural consequennce. The only exeption to this rule, would be a very special case, where $a(t)=t$, and k could non zero (like for example, the Milne universe). This will be discussed later. \\ For this metric, with $k=0$, the scale parameter should be $ a(t)=t^{\alpha} $, when $ \alpha $ is the index of the scale parameter: $ \alpha=\frac{1}{2} $
  for radiation, and $ \alpha=\frac{2}{3} $
  for dust (dark matter). With the homothetic vector filed of the metric, using (4) we get:
\begin{equation}
\chi^{\mu}=\frac{\chi}{4}\left(t,(1-\alpha)\overrightarrow{x}\right)
\end{equation}
 For this space we obtain that: $ \chi_{\mu;\nu}=\frac{\chi}{4}g_{\mu\nu} $ with no anti-symmetrical contribution to the tensor $ \chi_{\mu;\nu} $.
\section{Radiation Like Scalar Field}

Some approches in the context of Two Measure Theories [3] have shown that a scalar field can act as dust (dark matter). Here we will see that in our approach a scalar field can have the behavior of radiation. We choose a stress energy tensor of scalar filed $ \phi $, which doesn't have a trace $ T(\chi)=0 $, therefore we choose:

\begin{equation}
T_{(\chi)}^{\mu\nu}=(g^{\mu\alpha}g^{\nu\beta}-\frac{1}{4}g^{\mu\nu}g^{\alpha\beta})\phi_{,\alpha}\phi_{,\beta}
\end{equation}

From the Cosmological Principle the dependence of the scalar filed is only a function of time $ \phi(t) $. From (7) we get $ T_{\nu(\chi)}^{\mu}=\rho_{(\chi)}diag\left\{ 1,-1/3,-1/3,-1/3\right\}  $ when $ \rho(\chi)=(\frac{d}{dt}\phi)^{2} $. From conservation of $ T_{\left(\chi\right)}^{\mu\nu} $
 we get $ \rho(\chi)\sim a^{-4} $, like radiation. From variation of $ \phi $ we get a covariantly conserved current:

\begin{equation}
j_{(\chi)}^{\mu}=\left(\chi^{\mu;\nu}+\chi^{\nu;\mu}\right)\phi_{,\nu}-\frac{1}{2}\chi\phi^{,\mu}
\end{equation}

The conservation of the current $ j_{;\mu}^{\mu}(\chi)=0 $
  may be a condition for the vector field $ \chi^{\lambda} $. For the ansatz (4) of homothetic Killing vector we get $ j_{(\chi)}^{\mu}=0 $. This obviously satisfies the covariant conservation of (8). From the variation of the metric we get the effective stress-energy momentum tensor,which is the right hand side of the Einstein Tensor, then (5) reduse to the term:

\begin{equation}
T_{eff}^{\mu\nu}=-\frac{\chi}{2}T_{(\chi)}^{\mu\nu}-\frac{1}{2}\left(\chi^{\lambda}\nabla_{\lambda}\right)T_{(\chi)}^{\mu\nu}=-\frac{1}{2}(\chi^{\lambda}T_{(\chi)}^{\mu\nu})_{;\lambda}
\end{equation}

Two Friedman equation's are reduced to one, because of the tracelessness of the stress-energy tensor. Then we take the '00' compoent for the and $\rho(\chi)\sim a^{-4} $, then we get:

\begin{equation}
\rho_{eff}=-\frac{\chi}{2}(1+\frac{1}{2}t\frac{\partial}{\partial t})\rho_{(\chi)}=\chi\frac{\alpha-1}{2}\rho_{(\chi)}
\end{equation}

when Hubble Constant is defined as $ H=\frac{\dot{a}}{a} $
  (for a power law universe $ H=\alpha/t $). That means $ \rho_{eff} $ is proportional to $ \rho(\chi) $. We can get from that $ a\propto t^{1/2} $,because of the relation $ \rho_{eff}\propto H^{2} $, which cames from Einstein equation, for FRLW metric. The constans factor between is with the constant factor between $ \rho_{eff} $ and $ \rho(\chi) $ is $ -\frac{\chi}{4}$. We understand that the physical meaning of this theory is that a scalar filed could act like a radiation field, and could co-exist with the backround radiation in the early universe. The question we would like to search is do this model of scalar filed can predict any existence of dust or dark energy, if there is any potential for this field.

 \section{Electromagnetic Radiation}
  Another possibilty to explore is when $ T_{(\chi)}^{\mu\nu} $ will be like the standard Electromagnetic stress-energy mommentum tensor:

 \begin{equation}
  T_{(\chi)}^{\mu\nu}=F^{\mu\alpha}F_{\alpha}^{\nu}-\frac{1}{4}g^{\mu\nu}F_{\alpha\beta}F^{\alpha\beta}
 \end{equation}

  where $ F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu} $. From a variation of $ A_{\mu} $  we get some anti-symmetric tensor, which can be interpreted as a “dielectric” field tensor:

  \begin{equation}
D^{\mu\nu}=-\chi^{\mu;\beta}F_{\beta}^{\nu}+\chi^{\nu;\beta}F_{\beta}^{\mu}-\chi^{\beta;\nu}F_{\beta}^{\mu}+\chi^{\beta;\mu}F_{\beta}^{\nu}-\chi F^{\mu\nu}
 \end{equation}

  which has no sources $ (\nabla_{\mu}D^{\mu\nu}=0
   ) $. For the ansatz (4) of homothetic Killing vector we get $ D^{\mu\nu}=0 $
   , like the effective current of the radiation like scalar field. In this case we will get the same identity like (11), and the field in dynamical time will give the same predictions abuot the behavior of the expantion of the universe. Unlike the case of a scalar filed, which is automatically homogeneous and isotropic, now for the $ A_{\nu} $
    field, we obtain that in order to regain homogenity and isotropy, we have to consider a statistical average.

  More general theory, without any scalar fields, may include the traditional action for E.M field and the action from (3). So the total action will be:

\begin{equation}
  S=\int d^{4}x[\sqrt{-g}\chi_{\mu;\nu}T_{\left(\chi\right)}^{\mu\nu}+\frac{1}{16\pi G}\sqrt{-g}R-\frac{1}{4}\sqrt{-g}F_{\alpha\beta}F^{\alpha\beta}]
\end{equation}

  Without the last term $ F_{\alpha\beta}F^{\alpha\beta} $
    , Maxwell equations are not obtained, even for spaces where the ansatz (4) of homothetic Killing vector applies, although the energy momentum conservation of $ T_{(\chi)}^{\mu\nu} $
    is obtained. If the space does not satisfies the ansatz (4), then there will be a deviation from Maxwell equations. This deviation will take the following form:
  \begin{equation}
D^{\mu\nu}=-\chi^{\mu;\beta}F_{\beta}^{\nu}+\chi^{\nu;\beta}F_{\beta}^{\mu}-\chi^{\beta;\nu}F_{\beta}^{\mu}+\chi^{\beta;\mu}F_{\beta}^{\nu}-\chi F^{\mu\nu}-F^{\mu\nu}
 \end{equation}

    where $\nabla_{\mu}D^{\mu\nu}=0 $. If the space satisfies the ansatz (4) we get the regular Maxwell equations. When the space time does not satisfy the ansatz (4), it aquires non-trivial diaelectric properties, that can exist close to charged massive body or black holes for example, which will change the general behavior of electromanetic field in curved space-time. The stress energy tensor for the case of homothetic Killing vector is:
    \begin{equation}
  T_{eff}^{\mu\nu}=T_{(\chi)}^{\mu\nu}-\frac{1}{2}(\chi^{\lambda}T_{(\chi)}^{\mu\nu})_{;\lambda}
  \end{equation}
  For non-abelian Yang-Mills gauge theories the same effect takes place, where $ F_{\mu\nu}^{a}=\partial_{\mu}A_{\nu}^{a}-\partial_{\nu}A_{\mu}^{a}+gf^{abc}A_{\mu}^{b}A_{\nu}^{c} $ , we will get that similar results analogous to those considered in the abelian case.
\section{Situation for which $T_{eff}^{\mu\nu}=0$ but $ T_{\left(\chi\right)}^{\mu\nu}\neq0$}
Looking for example at (15), with respect to cosmological solutions, like FRWM we get $ \rho_{eff}=[1+\chi\frac{\alpha-1}{2}]\rho_{(\chi)}$. In a particular case, the effective gravitational energy is zero $T_{eff}^{\mu\nu}=0$, but it contains radiation which is represented by $ T_{\left(\chi\right)}^{\mu\nu} $ and $\chi=4$. In this case we get that k doesn't have to be zero, in order to have the non trivial homogenity $\frac{1}{2}\left(\chi^{\lambda}\nabla_{\lambda}\right)T_{(\chi)}^{\mu\nu}=nT_{(\chi)}^{\mu\nu} $ which we discused at (5). Taking the case of Milne Universe ($k=-1$ and $a(t)\sim t$), and taking the case of homothetic Killing vector which is of the form $\chi^{\mu}=c(t,0,0,0)$, we obtained that $\rho_{(\chi)}\sim a^{-4}$, and because it's a flat universe, and  $\rho_{(\chi)}\sim t^{-4}$, which is a boundery between inflationary behavior ($\alpha>1$) and non-inflationary ($1>\alpha>0$) behavior.
\section{Discussion, Conclusions and Prospects}

   We have found that a generalization of TMT gives rise to two non trivial energy momentum tensors which are related but not equal in general. Because of the vector field variation, it enforces its covariant conservation of $ T_{\left(\chi\right)}^{\mu\nu} $
     in addition to $ T_{eff}^{\mu\nu} $
    . We have studied this for the case of the vector field as homothetic Killing vector, and using $ T_{\left(\chi\right)}^{\mu\nu} $
     of radiation or radiation like scalar filed Cosmological solutions. In the radiation era of the universe, the $ T_{eff}^{\mu\nu} $
     stress energy tensor could co-exist with the original radiation, which is the source of the Backround Cosmic Radiation. The solution, which could be a model for the early universe, requires no spacial curvature for the universe, because those solutions require homogenity of the dependence on time of the scale factor $ a(t)=t^{\alpha} $
     in order for the $ T_{\left(\chi\right)}^{\mu\nu} $
     and $ T_{eff}^{\mu\nu} $
     to be proportional to each other. This mathematical explanation singles out zero spacial curvature. So it may be an alternative explanation for negligible spacial curvature instead of inflation [4]. In any case, inflation still takes care of other problems, like the horizon problem.

   These kind of models could be used in the context of inlation. It will be interesting to add a potential to this scalar filed, and see if one can reproduce an inflationary behavior, also a quantum behavior of this field, and its affect on the first moments of the universe as an inflationary model.

   When $ T_{\left(\chi\right)}^{\mu\nu} $
     is a standard stress energy tensor of gauge filed, we obtain similar conclusions. Also this could lead to a modified electromagnetism and similar situations for other gauge theories in curved space time. We would like to find the spherically symmetric charged solutions like charged black hole.

     \section{Appendix - identities}
     $$\frac{\partial g^{\alpha\beta}}{\partial g_{\mu\nu}}=-\frac{1}{2}(g^{\alpha\mu}g^{\beta\nu}+g^{\alpha\nu}g^{\beta\mu})$$\\$$\frac{\partial\Gamma_{\lambda\sigma}^{\tau}}{\partial g_{\mu\nu}}=-\frac{1}{2}(g^{\mu\tau}\Gamma_{\lambda\sigma}^{\nu}+g^{\nu\tau}\Gamma_{\lambda\sigma}^{\mu})$$\\$$\frac{\partial\Gamma_{\lambda\alpha}^{\tau}}{\partial g_{\mu\nu,\sigma}}=\frac{1}{4}\left[g^{\mu\tau}\left(\delta_{\alpha}^{\nu}\delta_{\lambda}^{\sigma}+\delta_{\lambda}^{\nu}\delta_{\alpha}^{\sigma}\right)+g^{\mu\nu}\left(\delta_{\alpha}^{\mu}\delta_{\lambda}^{\sigma}+\delta_{\lambda}^{\mu}\delta_{\alpha}^{\sigma}\right)-g^{\tau\sigma}\left(\delta_{\alpha}^{\mu}\delta_{\lambda}^{\nu}+\delta_{\lambda}^{\mu}\delta_{\alpha}^{\nu}\right)\right]$$
     \\$$T_{eff}^{\alpha\beta}=\frac{1}{\sqrt{-g}}\frac{\partial\left(\sqrt{-g}\chi_{\mu;\nu}T_{\left(\chi\right)}^{\mu\nu}\right)}{\partial g_{\alpha\beta}}-\frac{1}{\sqrt{-g}}\frac{\partial}{\partial x^{\sigma}}\frac{\partial\left(\sqrt{-g}\chi_{\mu;\nu}T_{\left(\chi\right)}^{\mu\nu}\right)}{\partial g_{\alpha\beta,\sigma}}$$


\section{Introduction}

Nd-Fe-B-based nanocomposite permanent magnets, which consist of exchange-coupled nanocrystalline hard ($\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}$) and soft ($\alpha$-$\mathrm{Fe}$ or $\mathrm{Fe}_3\mathrm{B}$) magnetic phases, are of potential interest for electronic devices due to their preeminent magnetic properties such as high remanence and magnetic energy product \cite{gutfleisch2011,liu2009}. The major challenge remains the understanding of how the details of the microstructure (e.g., average particle size and shape, volume fraction of soft phase, texture, interfacial chemistry) correlate with their magnetic properties. In order to tackle this issue a multiscale characterization approach is adopted, which comprises a suite of both experimental and theoretical state-of-the-art methods such as high-resolution electron microscopy, electron backscattering diffraction, three-dimensional atom-probe analysis, Lorentz and Kerr microscopy, or atomistic and continuum micromagnetic simulations.

Recent investigations by Liu \textit{et al.} \cite{liu2013} and Sepehri-Amin \textit{et al.} \cite{SepehriAmin2013} demonstrate that the properties of the interface regions between the $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}$ grains decisively determine the coercivity of the sample. The grain-boundary layers (and triple-junctions between the grains) have a thickness between about $1-15 \, \mathrm{nm}$ and can be both in a crystalline or amorphous state. Moreover, as far as their magnetism is concerned, the intergranular regions are characterized by different magnetic interactions (exchange, magnetocrystalline anisotropy, saturation magnetization) as compared to the $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}$ crystallites and, hence, they represent potential sources for the nucleation of inhomogeneous spin textures during magnetization reversal. Indeed, the micromagnetic simulation results reported in \cite{SepehriAmin2013} suggest that the existence of a thin ($< 5 \, \mathrm{nm}$) ferromagnetic grain-boundary phase with reduced magnetocrystalline anisotropy, exchange-stiffness constant, and saturation magnetization causes the magnetization reversal to occur from the soft intergranular phase into the hard $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}$ phase at a field of $-2.5 \, \mathrm{T}$. When the intergranular phase is nonmagnetic, then the nucleation of reversed domains starts from the triple junctions of the $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}$ grains at a higher field of $-3.2 \, \mathrm{T}$.

In this context it also worth noting that first-principles density functional theory calculations on an exchange-spring multilayer system \cite{toga2016} predict a dependency of the exchange coupling on the crystallographic orientation at the interface between $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}$ and $\alpha$-$\mathrm{Fe}$; specifically, ferromagnetic coupling is predicted for the $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}(001)/$$\alpha$-$\mathrm{Fe}(001)$ interface model, whereas antiferromagnetic interactions are obtained for $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}(100)/$$\alpha$-$\mathrm{Fe}(110)$. If this prediction were true, then it may negatively influence the magnetic properties of this class of materials (e.g., the maximum energy product). Indeed, a recent experimental study using ferromagnetic resonance and Kerr microscopy \cite{ogawa2015} reports on the predicted negative exchange coupling.

The above discussed examples ultimately demonstrate that the magnetic microstructure of nanocrystalline Nd-Fe-B-based magnets is characterized by inhomogeneous magnetization structures and that interfacial regions are a major cause for the nanoscale spin disorder. In addition to the grain boundaries, there exist, however, other sources of spin disorder in such materials: ultrafine-grained textured nanocomposites are produced from melt-spun ribbons via hot compaction \cite{gutfleisch2000,kirchner2004,liu2013,SepehriAmin2013}; this process may introduce crystal defects which locally act as nucleation centers for nonuniform magnetization textures. Furthermore, one has to invoke a magnetization inhomogeneity which is due to the nonideal alignment (dispersion) of the crystallographic $c$-axes (of the $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}$ grains) along the pressing direction during hot deformation; the spins have to undergo rotations in order to accommodate to the changes in the easy-axis magnetization direction from grain to grain. Last but not least, there is the magnetic shape anisotropy of the usually platelet-shaped $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}$ particles, which may result in a small spin canting towards the plane perpendicular to the $c$-axis. It is certainly true that the magnetic anisotropy field of the $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}$ phase (about $8 \, \mathrm{T}$ at $300 \, \mathrm{K}$ \cite{woodcock2012}) is much larger than any shape-anisotropy field (assuming, e.g., $0.5 \, \mathrm{T}$ for strongly anisotropic grains), but nevertheless weak spin canting ($\tan^{-1}(0.5/8) \cong 3.6^{\circ}$) might be produced by the competition between shape and magnetocrystalline anisotropy.

In order to scrutinize the above-sketched issue, we have carried out a comparative study of the magnetic microstructure of textured (hot-deformed) and isotropic nanocrystalline $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}/\alpha$-$\mathrm{Fe}$ by means of magnetic small-angle neutron scattering (SANS). Specifically, the central aim of our investigation is to detect and quantify the presumed nanoscale spin disorder, which is commonly only indirectly inferred by combining results from electron microscopy, magnetization, and micromagnetic simulations.

The SANS technique (see Ref.~\cite{michels2014review} for a review) provides information on variations of both the magnitude and orientation of the magnetization on a nanometer length scale ($\sim 1-300$~nm). SANS is extremely sensitive to long-wavelength magnetization fluctuations and it has only recently been employed for characterizing Nd-Fe-B-based permanent magnets: for example, the field dependence of characteristic magnetic length scales during the magnetization-reversal process in isotropic Nd-Fe-B-based nanocomposites \cite{bickapl2013} and in isotropic sintered Nd-Fe-B \cite{perigo2015} was studied, the exchange-stiffness constant has been determined \cite{bickapl2013no2}, the observation of the so-called spike anisotropy in the magnetic SANS cross section has been explained with the formation of flux-closure patterns \cite{perigo2014}, magnetic multiple scattering has been detected \cite{ono2016}, textured Nd-Fe-B has been investigated \cite{perigo2016no1}, and the effect of grain-boundary diffusion on the magnetization-reversal process of isotropic \cite{perigo2016no2} and hot-deformed textured \cite{yano2012,yano2014,saito2015} nanocrystalline Nd-Fe-B magnets has been studied.

\section{Experimental}

Two $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}/\alpha$-$\mathrm{Fe}$ nanocomposites containing, respectively, $5 \, \mathrm{wt\%}$ of Fe were investigated in this study. Both samples were prepared by means of the melt-spinning technique. One sample was subsequently hot-deformed in order to obtain a textured magnet. For this purpose, the melt-spun ribbons were crushed into powders of a few hundred micrometers in size and then sintered at $973 \, \mathrm{K}$ under a pressure of $400 \, \mathrm{MPa}$. The sintered bulk was hot-deformed with a height reduction of about $75 \, \%$ to develop the [001] texture of the $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}$ phase along the pressing direction \cite{kirchner2004,SepehriAmin2013}. This results in the formation of platelet-shaped $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}$ grains with an average thickness of $\sim 110 \, \mathrm{nm}$ and an average diameter of $\sim 140 \, \mathrm{nm}$. The $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}$ platelets are stacked along the nominal $c$-axis, which we define as the [001] direction, with some degrees of misorientation. The isotropic sample had an average grain size of about $20 \, \mathrm{nm}$. We have also investigated composites with $0 \, \mathrm{wt\%}$ and $10 \, \mathrm{wt\%}$ of Fe, which, as far as the neutron results are concerned, show qualitatively the same behavior as the $5 \, \mathrm{wt\%}$ sample. For further details, see Refs.~\cite{yano2012,yano2014,saito2015}.

The neutron experiment has been carried out at $300 \, \mathrm{K}$ at the instrument SANS-I at the Paul Scherrer Institute, Switzerland, using unpolarized neutrons with a mean wavelength of $\lambda = 4.5 \, \mathrm{\AA}$ and $\Delta \lambda / \lambda = 10 \, \mathrm{\%}$ (FWHM) \cite{kohlbrecher2000,niketic2015}. The external magnetic field $\mathbf{H}_0$ (provided by a cryomagnet; $\mu_0 H_{\mathrm{max}} = 9.5 \, \mathrm{T}$) was applied perpendicular and parallel to the wave vector $\mathbf{k}_0$ of the incoming neutron beam (compare Fig.~\ref{fig1}); this corresponds to the situation that $\mathbf{H}_0$ is parallel ($\mathbf{k}_0 \perp \mathbf{H}_0$) and perpendicular ($\mathbf{k}_0 \parallel \mathbf{H}_0$) to the nominal $c$-axis (pressing direction) of the textured sample. Neutron data were corrected for background scattering (empty sample holder), transmission, and detector efficiency using the GRASP software package. The measured transmission was larger than $90 \, \mathrm{\%}$ for both samples at all fields investigated. Further sample characterization was done by means of vibrating sample magnetometry, transmission electron microscopy, and synchrotron x-ray diffraction (at beamline BL04-MSPD at the Alba synchrotron, Barcelona, Spain \cite{Fauth2015}).

\begin{figure}[htb!]
\includegraphics[width=1.0\columnwidth]{fig1.eps}
\caption{\label{fig1} Sketch of the perpendicular (a) and parallel (b) scattering geometry, which, respectively, have the applied magnetic field $\mathbf{H}_0$ perpendicular and parallel to the wave vector $\mathbf{k}_0$ of the incident neutron beam; $q = \left|\mathbf{q}\right| = 4\pi \lambda^{-1} \sin\psi$, where $2\psi$ denotes the scattering angle and $\lambda$ is the mean neutron wavelength. Note that $\mathbf{H}_0 \parallel \mathbf{e}_z$ in \textit{both} geometries and that $\mathbf{q} \cong (0, q_y, q_z) = q (0, \sin\theta, \cos\theta)$ for $\mathbf{k}_0 \perp \mathbf{H}_0$ (a) and $\mathbf{q} \cong (q_x, q_y, 0) = q (\cos\theta, \sin\theta, 0)$ for $\mathbf{k}_0 \parallel \mathbf{H}_0$ (b). The pressing direction is horizontal, which is along $\mathbf{e}_z$ in (a) and along $\mathbf{e}_x$ in (b).}
\end{figure}

\section{Unpolarized SANS cross sections and correlation function}

The elastic unpolarized SANS cross section $d \Sigma / d \Omega$ at momentum-transfer vector $\mathbf{q}$ takes on different forms depending on the relative orientation between the wave vector $\mathbf{k}_0$ of the incident neutron beam and the externally applied magnetic field $\mathbf{H}_0$ \cite{michels2014review}; for the perpendicular geometry  ($\mathbf{k}_0 \perp \mathbf{H}_0$), we obtain
\begin{eqnarray}
\label{sigmatotperp}
\frac{d \Sigma}{d \Omega}(\mathbf{q}) = \frac{8 \pi^3}{V} \left( |\widetilde{N}|^2 + b_H^2 |\widetilde{M}_x|^2 + b_H^2 |\widetilde{M}_y|^2 \cos^2\theta \right. \nonumber \\ \left. + b_H^2 |\widetilde{M}_z|^2 \sin^2\theta - b_H^2 (\widetilde{M}_y \widetilde{M}_z^{\ast} + \widetilde{M}_y^{\ast} \widetilde{M}_z) \sin\theta \cos\theta \right) , \nonumber \\
\end{eqnarray}
whereas for the parallel case ($\mathbf{k}_0 \parallel \mathbf{H}_0$)
\begin{eqnarray}
\label{sigmatotpara}
\frac{d \Sigma}{d \Omega}(\mathbf{q}) = \frac{8 \pi^3}{V} \left( |\widetilde{N}|^2 + b_H^2 |\widetilde{M}_x|^2 \sin^2\theta + b_H^2 |\widetilde{M}_y|^2 \cos^2\theta \right. \nonumber \\ \left. + b_H^2 |\widetilde{M}_z|^2 - b_H^2 (\widetilde{M}_x \widetilde{M}_y^{\ast} + \widetilde{M}_x^{\ast} \widetilde{M}_y) \sin\theta \cos\theta \right) ; \nonumber \\
\end{eqnarray}
$V$ denotes the scattering volume, $b_H = 2.91 \times 10^{8} \, \mathrm{A^{-1} m^{-1}}$, $\widetilde{N}(\mathbf{q})$ is the nuclear scattering amplitude, and $\widetilde{\mathbf{M}}(\mathbf{q}) = \{ \widetilde{M}_x(\mathbf{q}), \widetilde{M}_y(\mathbf{q}), \widetilde{M}_z(\mathbf{q}) \}$ represents the Fourier transform of the magnetization $\mathbf{M}(\mathbf{r}) = \{ M_x(\mathbf{r}), M_y(\mathbf{r}), M_z(\mathbf{r}) \}$; $c^{*}$ is a quantity complex-conjugated to $c$. We would like to emphasize that the magnetization vector field of a bulk ferromagnet is a function of the position $\mathbf{r} = \{ x, y, z \}$ inside the material, i.e., $\mathbf{M} = \mathbf{M}(x, y, z)$, and that, consequently, $\widetilde{\mathbf{M}} = \widetilde{\mathbf{M}}(q_x, q_y, q_z)$. However, the Fourier components which appear in the above SANS cross sections represent projections into the $q_y$-$q_z$-plane for $\mathbf{k}_0 \perp \mathbf{H}_0$ ($q_x \cong 0$) and into the $q_x$-$q_y$-plane for $\mathbf{k}_0 \parallel \mathbf{H}_0$ ($q_z \cong 0$) (compare Fig.~\ref{fig1}). In polar coordinates, the $\widetilde{M}_{x,y,z}$ then depend (in addition to the applied field and the magnetic interactions) on both the magnitude $q$ and the orientation $\theta$ of the scattering vector $\mathbf{q}$ \cite{erokhin2015}.

In our neutron data analysis below, we subtract the respective SANS signal at the largest available field of $9.5 \,  \mathrm{T}$ (approach-to-saturation regime, compare Fig.~\ref{fig2}) from the measured data at lower fields. This subtraction procedure eliminates the nuclear SANS contribution ($\propto |\widetilde{N}|^2$), which is field independent, and it yields the so-called spin-misalignment SANS cross section $d \Sigma_M / d \Omega$, which we display here for simplicity only for the parallel scattering geometry:
\begin{eqnarray}
\label{sigmasmpara}
\frac{d \Sigma_M}{d \Omega} = \frac{8 \pi^3}{V} \, b_H^2 \left( \Delta |\widetilde{M}_x|^2 \sin^2\theta + \Delta |\widetilde{M}_y|^2 \cos^2\theta \right. \nonumber \\ \left. + \Delta |\widetilde{M}_z|^2 + \Delta CT \sin\theta \cos\theta \right) ,
\end{eqnarray}
where $\Delta |\widetilde{M}_x|^2 := |\widetilde{M}_x|^2(H) - |\widetilde{M}_x|^2(9.5 \,  \mathrm{T})$ (and so on for the other Fourier coefficients) represents the difference between the value of $|\widetilde{M}_x|^2$ at the actual field $H$ and the measurement at $9.5 \,  \mathrm{T}$ [$CT := - (\widetilde{M}_x \widetilde{M}_y^{\ast} + \widetilde{M}_x^{\ast} \widetilde{M}_y$)]. If it would be possible to fully saturate the sample (i.e., $\mathbf{M}(\mathbf{r}) = \{ 0, 0, M_z = M_s(\mathbf{r}) \}$) and if one restricts the considerations (subtraction procedure) to the approach-to-saturation regime, where the field dependence of the longitudinal Fourier component can be neglected (i.e., $|\widetilde{M}_z|^2(H) - |\widetilde{M}_s|^2 \rightarrow 0$), then $d \Sigma_M / d \Omega$ (for $\mathbf{k}_0 \parallel \mathbf{H}_0$) reduces to
\begin{eqnarray}
\label{sigmasmparainfinity}
\frac{d \Sigma_M}{d \Omega} = \frac{8 \pi^3}{V} \, b_H^2 \left( |\widetilde{M}_x|^2 \sin^2\theta + |\widetilde{M}_y|^2 \cos^2\theta \right. \nonumber \\ \left. + CT \sin\theta \cos\theta \right) ,
\end{eqnarray}
and likewise for the $\mathbf{k}_0 \perp \mathbf{H}_0$ geometry.

Furthermore, it is decisive for the later discussion to note that for a ferromagnet with a statistically-isotropic microstructure the parallel total (nuclear and magnetic) $d \Sigma / d \Omega$ and $d \Sigma_M / d \Omega$ [Eqs.~(\ref{sigmatotpara}) and (\ref{sigmasmpara})] are generally \textit{isotropic}, i.e., $\theta$-independent (see, e.g., Fig.~21 in Ref.~\cite{michels2014review}, Fig.~4 in Ref.~\cite{michels2014jmmm}, or Figs.~\ref{fig5}(b), \ref{fig7}(c), and \ref{fig7}(d) below). In other words, although the individual contributions to the parallel SANS cross section are highly anisotropic (e.g., $|\widetilde{M}_x|^2 \sin^2\theta$), their corresponding sums in Eqs.~(\ref{sigmatotpara}) and (\ref{sigmasmpara}) are isotropic for a statistically-isotropic ferromagnet; this is not true for the perpendicular geometry [Eq.~(\ref{sigmatotperp})], which generally exhibits a pronounced angular anisotropy.

The (normalized) correlation function $c(r)$ of the spin misalignment can be computed from azimuthally-averaged data via \cite{mettus2015}
\begin{equation}
\label{corrfunc}
c(r) = \frac{\int_0^{\infty} \frac{d \Sigma_M}{d \Omega}(q) J_0(q r) \, q \, dq}{\int_0^{\infty} \frac{d \Sigma_M}{d \Omega}(q) \, q \, dq} ,
\end{equation}
where $J_0(qr)$ denotes the zeroth-order Bessel function. Analysis of $c(r)$ provides information on the characteristic magnetic length scales \cite{bickapl2013,perigo2015,perigo2016no2}.

\begin{figure*}[tb!]
\includegraphics[width=1.65\columnwidth]{fig2.eps}
\caption{\label{fig2} Room-temperature magnetization curves of textured (a) and isotropic (b) $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}/\alpha$-$\mathrm{Fe}$ ($5 \, \mathrm{wt\%}$ Fe). Measurements have been carried out for the magnetic field applied parallel and perpendicular to the texture axis (pressing direction) in (a), and for two different in-plane directions in (b) (``in-plane~2'' direction is rotated by $90^{\circ}$ with respect to ``in-plane~1'' direction). Magnetization data (on the rectangular-shaped samples) have been corrected for demagnetizing effects using the magnetometric demagnetizing factor \cite{aharoni98}.}
\end{figure*}

\begin{figure*}[tb!]
\centering
\includegraphics[width=1.90\columnwidth]{fig3.eps}
\caption{\label{fig3} Bright-field transmission electron microscopy images of the textured [(a) and (b)] and isotropic (c) $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}/\alpha$-$\mathrm{Fe}$ nanocomposites ($5 \, \mathrm{wt\%}$ Fe). The average sizes of the anisotropic grains of the textured sample have been estimated [from (a) and (b)] as, respectively, $\sim 110 \, \mathrm{nm}$ (parallel to the pressing direction ``$p$'') and $\sim 140 \, \mathrm{nm}$ (perpendicular to the pressing direction ``$p$''), while the average grain diameter of the isotropic sample has been found to be $\sim 20 \, \mathrm{nm}$.}
\end{figure*}

\begin{figure*}[htb!]
\centering
\includegraphics[width=2.0\columnwidth]{fig4.eps}
\caption{\label{fig4} Synchrotron x-ray diffraction data of the isotropic (a) and textured (b) $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}/\alpha$-$\mathrm{Fe}$ nanocomposites ($5 \, \mathrm{wt\%}$ Fe). The pressing direction of the hot-deformed sample is horizontal (same scattering geometry as in the neutron experiment, compare Fig.~\ref{fig1}). (left images) Integrated intensity as a function of azimuthal angle $\theta$ and scattering angle $2 \psi$; (right images) corresponding Debye-Scherrer diffraction rings. Radial integration of synchrotron data has been performed with the Fit2D software \cite{fit2d}.}
\end{figure*}

\section{Results and discussion}

Figure~\ref{fig2} displays the room-temperature magnetization curves of textured [Fig.~\ref{fig2}(a)] and isotropic [Fig.~\ref{fig2}(b)] $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}/\alpha$-$\mathrm{Fe}$. The coercive fields are $\mu_0 H_c = 0.57 \, \mathrm{T}$ (textured) and $\mu_0 H_c = 0.61 \, \mathrm{T}$ (isotropic). The saturation polarization was estimated by extrapolating the data to infinite field: we find $J_s = 1.57 \, \mathrm{T}$ (textured) and $J_s = 1.53 \, \mathrm{T}$ (isotropic) with ensuing remanence-to-saturation ratios of about $0.64$ (textured easy), $0.36$ (textured hard), and $0.53$ (isotropic). Consistent with the magnetization data, the transmission electron microscopy images of the textured sample [Fig.~\ref{fig3}(a) and \ref{fig3}(b)] reveal a weakly anisotropic microstructure, while the isotropic sample [Fig.~\ref{fig3}(c)] exhibits equiaxed grains.

X-ray diffraction measurements carried out (in transmission mode) at the Alba synchrotron (Fig.~\ref{fig4}) unambiguously prove the presence of a weak texture along the (horizontal) pressing direction; specifically, diffraction peaks of the type $(0 0 l)$ do present two maxima around $\theta = 0^{\circ}$ and $\theta = 180^{\circ}$ in the Debye Scherrer rings. Additionally, we find evidence for the presence of texture along other crystallographic directions; peaks of the type $(h k 0)$ exhibit two maxima around $90^{\circ}$ and $270^{\circ}$. This latter observation will be of relevance when discussing the results of the magnetic neutron data analysis (see below).

Figure~\ref{fig5} depicts (for $\mathbf{k}_0 \parallel \mathbf{H}_0$) the two-dimensional unpolarized total scattering cross sections $d \Sigma / d \Omega$ of the textured and isotropic Nd-Fe-B-based nanocomposites at selected applied magnetic fields ($9.5 \,  \mathrm{T}$, remanence, coercive field). The isotropic sample [Fig.~\ref{fig5}(b)] exhibits an isotropic scattering pattern at all fields investigated, whereas the textured sample [Fig.~\ref{fig5}(a)] shows anisotropic scattering with an elongation along the horizontal direction. The corresponding (over $2\pi$) azimuthally-averaged data sets are displayed in Fig.~\ref{fig6}; between the coercive field and the largest available field of $9.5 \,  \mathrm{T}$, the cross section of the isotropic sample changes (roughly) by about an order of magnitude at the smallest momentum-transfers $q$ (and about half an order of magnitude for the textured sample).

\begin{figure*}[tb!]
\includegraphics[width=1.90\columnwidth]{fig5.eps}
\caption{\label{fig5} Color-coded two-dimensional intensity maps of the total unpolarized $d \Sigma / d \Omega$ in the plane perpendicular to the incoming neutron beam at selected applied magnetic fields (see insets) (logarithmic color scale) ($\mathbf{k}_0 \parallel \mathbf{H}_0$). $d \Sigma / d \Omega$ of the textured (a) and isotropic (b) $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}/\alpha$-$\mathrm{Fe}$ nanocomposite. $\mathbf{H}_0$ is normal to the detector plane.}
\end{figure*}
\begin{figure*}[tb!]
\includegraphics[width=2.0\columnwidth]{fig6.eps}
\caption{\label{fig6} Azimuthally-averaged total unpolarized SANS cross sections $d \Sigma / d \Omega$ at selected applied magnetic fields (see insets) (log-log scale) ($\mathbf{k}_0 \parallel \mathbf{H}_0$). $d \Sigma / d \Omega$ of the textured (a) and isotropic (b) $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}/\alpha$-$\mathrm{Fe}$ nanocomposite. (c) Applied-field dependence of the power-law exponent $n$ in $d \Sigma_M / d \Omega = K/q^n$ for the textured and isotropic $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}/\alpha$-$\mathrm{Fe}$ nanocomposite. $d \Sigma_M / d \Omega$ has been obtained by subtracting, respectively, the total $d \Sigma / d \Omega$ at $9.5 \, \mathrm{T}$; the fits were restricted to the interval $0.4 \, \mathrm{nm}^{-1} \lesssim q \lesssim 0.6 \, \mathrm{nm}^{-1}$. Dotted horizontal line ($n = 4$) corresponds to scattering due to sharp interfaces (Porod) or to exponentially correlated magnetization fluctuations.}
\end{figure*}

While the textured nanocomposite reveals a power-law type scattering over most of the $q$-range, the isotropic sample exhibits a more structured $d \Sigma / d \Omega$ with significant curvature at lower and medium $q$. This difference in $d \Sigma / d \Omega$ is most likely related to the difference in the average grain sizes and the ensuing magnetization fluctuations on a nanometer length scale: the isotropic sample has an average grain size of $\sim 20 \, \mathrm{nm}$, while the textured Nd-Fe-B possesses a larger particle size of the order of $100 \,\mathrm{nm}$ (compare the TEM images in Fig.~\ref{fig3}). We also note that the $d \Sigma / d \Omega$ of both samples (data not shown) as well as the spin-misalignment SANS cross section $d \Sigma_M / d \Omega$ [Fig.~\ref{fig6}(c)] are characterized by power-law exponents $n$ that are larger than $4$. This is in agreement with the notion of spin-misalignment scattering, i.e., scattering due to canted spins with a characteristic magnetic-field-dependent wavelength \cite{michels2014review}. It is also quite obvious from this observation that the corresponding magnetization fluctuations in real space are not exponentially correlated (see Fig.~\ref{fig8} below).

As discussed previously, for $\mathbf{k}_0 \parallel \mathbf{H}_0$, any anisotropy of $d \Sigma / d \Omega$ (or of $d \Sigma_M / d \Omega$) is indicative of an anisotropic microstructure. At magnetic \textit{saturation}, the \textit{total} SANS signal arises from nanoscale spatial fluctuations in the nuclear density and in the saturation magnetization $M_s(\mathbf{r})$, presumably at internal $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}/\alpha$-$\mathrm{Fe}$ interfaces. The nuclear scattering-length density contrast between the $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}$ phase and the $\alpha$-Fe phase amounts to $\Delta \rho_{\mathrm{nuc}} \cong 1.63 \times 10^{14} \, \mathrm{m}^{-2}$, whereas -- at saturation -- the magnetic contrast can be estimated as $\Delta \rho_{\mathrm{mag}} = b_H \Delta M \cong 1.37 \times 10^{14} \, \mathrm{m}^{-2}$, where $\Delta M$ denotes the difference in saturation magnetization between $\alpha$-Fe ($J_s = 2.2 \, \mathrm{T}$) and $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}$ ($J_s = 1.61 \, \mathrm{T}$). By assuming that the elements of the microstructure which give rise to nuclear scattering $|\widetilde{N}|^2$ are identical to those which give rise to longitudinal magnetic scattering $b_H^2 |\widetilde{M}_z|^2$, one finds for a saturated sample that the ratio of nuclear to magnetic SANS equals $|\widetilde{N}|^2/(b_H^2 |\widetilde{M}_z|^2) \cong 1.42$. With reference to the electron-microscopy results (Fig.~\ref{fig3}), which reveal a (weakly) anisotropic grain shape (aspect ratio $\sim 1.3$), it is then obvious that a (weakly) horizontally-elongated SANS pattern can already be generated \textit{at saturation} by the combined action of the nuclear and longitudinal magnetic form factors.

\begin{figure}
\includegraphics[width=1.0\columnwidth]{fig7.eps}
\caption{\label{fig7} Selected results for the spin-misalignment SANS cross section $d \Sigma_M / d \Omega$ of the textured and isotropic $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}/\alpha$-$\mathrm{Fe}$ nanocomposite for $\mathbf{k}_0 \parallel \mathbf{H}_0$ [(a)$-$(d)] and for $\mathbf{k}_0 \perp \mathbf{H}_0$ [(e)$-$(f)] (logarithmic color scale). The respective data set at the maximum applied field of $9.5 \, \mathrm{T}$ has been subtracted. In (a)$-$(d), $\mathbf{H}_0$ is normal to the detector plane, whereas in (e)$-$(f) $\mathbf{H}_0$ is horizontal in the plane.}
\end{figure}

Subtracting the total $d \Sigma / d \Omega$ at $9.5 \, \mathrm{T}$ from the total $d \Sigma / d \Omega$ at lower fields, we obtain the spin-misalignment SANS cross section $d \Sigma_M / d \Omega$ [Eq.~(\ref{sigmasmpara})], which is free of nuclear SANS. The results for $d \Sigma_M / d \Omega$ for the textured nanocomposite [Fig.~\ref{fig7}(a) and \ref{fig7}(b)] still reveal an angular anisotropy with maxima parallel and antiparallel to the horizontal texture axis. Inspection of Eq.~(\ref{sigmasmpara}) then suggests that this observation may be due to (i) spin components which are directed along the $\pm \mathbf{e}_y$-direction [cf.\ the term $\Delta |\widetilde{M}_y|^2 \cos^2\theta$ in Eq.~(\ref{sigmasmpara})] and/or due to (ii) the particle form factor anisotropy (cf.\ terms $\propto \Delta |\widetilde{M}_z|^2$). However, measurements in the $\mathbf{k}_0 \perp \mathbf{H}_0$ geometry [compare Fig.~\ref{fig1}(a)] suggest that longitudinal magnetization fluctuations play only a minor role: if the $d \Sigma_M / d \Omega$ (for $\mathbf{k}_0 \perp \mathbf{H}_0$) were dominated by $\Delta |\widetilde{M}_z|^2$, a $\sin^2\theta$-type anisotropy with intensity maxima along the vertical direction would result [compare Eq.~(\ref{sigmatotperp})]. This is, however, not visible in the experimental data [Fig.~\ref{fig7}(e) and \ref{fig7}(f)], which exhibit a horizontal elongation [cf.\ the term $|\widetilde{M}_y|^2 \cos^2\theta$ in Eq.~(\ref{sigmatotperp})]. In other words, the anisotropy of the scattering pattern for $\mathbf{k}_0 \parallel \mathbf{H}_0$ [Fig.~\ref{fig7}(a) and \ref{fig7}(b)] is due to an anisotropy in the magnetic microstructure, not to the form-factor anisotropy $\Delta |\widetilde{M}_z|^2$ of the particles.

We emphasize that, although the mean magnetization in the remanent state is directed along the $c$-axis [$\mathbf{e}_z$-direction in Fig.~\ref{fig1}(a) and $\mathbf{e}_x$-direction in Fig.~\ref{fig1}(b)], the magnetic neutron scattering cross section is in both geometries dominated by the respective $|\widetilde{M}_y|^2 \cos^2\theta$ term, which (in real space) is related to small misaligned spin components varying along the $\pm \mathbf{e}_y$-direction \cite{michels03epl}. This anisotropy in the spin microstructure may be related to the finding of a crystallographic texture: as shown in Fig.~\ref{fig4}, diffraction peaks of the type $(h k 0)$ exhibit two maxima along the vertical direction ($\theta = 90^{\circ}$ and $\theta = 270^{\circ}$). The investigation of the relation between this crystallographic texture and the spin texture is of interest in its own right and beyond the scope of this paper. However, we would like to emphasize that recent electron-microscopy and three-dimensional atom-probe tomography work by Liu \textit{et al.} \cite{liu2014} also reports anisotropic properties of the grain-boundary phase in hot-deformed nanocrystalline Nd-Fe-B magnets; namely, these authors found that the concentration of rare-earth elements is higher for intergranular phases parallel to the flat surface of the platelet-shaped $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}$ grains as compared to intergranular phases along the short side of the platelets.

The characteristic size of the spin inhomogeneities in the remanent state along the vertical and horizontal direction has been estimated by computing [using Eq.~(\ref{corrfunc})] the correlation function $c(r)$ of the spin misalignment (Fig.~\ref{fig8}). The $\exp(-1)$-lengths are $l_C \cong 53 \, \mathrm{nm}$ along the vertical direction, and $l_C \cong 42 \, \mathrm{nm}$ along the horizontal direction; $l_C \cong 28 \, \mathrm{nm}$ for the isotropic nanocomposite. Note that taking the $\exp(-1)$-lengths does not imply that the correlations decay exponentially. For the textured specimen, both $l_C$ values are smaller than the average particle size, which suggests the existence of intraparticle spin disorder, whereas $l_C \cong D $ for the isotropic sample. Compatible with \cite{liu2014}, these results indicate that the microscopic nature of the microstructural defects (e.g., the $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}/\alpha$-$\mathrm{Fe}$ interfaces) along these two directions are different (as is manifest by the different correlation lengths). In this respect field-dependent SANS measurements are helpful, since they allow one to determine the field evolution of $l_C$, from which the size of the defect (causing the spin perturbation) and the exchange correlation length may be obtained \cite{bickapl2013,perigo2015}.

\begin{figure}[tb!]
\includegraphics[width=0.80\columnwidth]{fig8.eps}
\caption{\label{fig8} Normalized correlation function $c(r)$ of the spin misalignment [Eq.~(\ref{corrfunc})] for the textured and isotropic $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}/\alpha$-$\mathrm{Fe}$ nanocomposite in the remanent state. $c(r)$ of the textured sample has been computed using $d \Sigma_M / d \Omega$ averaged along the vertical and horizontal directions ($\pm 7.5^{\circ}$ sector averages) as well as using the full circular ($2\pi$) average of $d \Sigma_M / d \Omega$; the $c(r)$ of the isotropic sample was computed using the corresponding $2\pi$-averaged $d \Sigma_M / d \Omega$ (see inset). Solid horizontal line: $C(r) = \exp(-1)$. The physically relevant information content of $c(r)$ is restricted to the interval $[r_{\mathrm{min}}, r_{\mathrm{max}}]$ with approximately $r_{\mathrm{min}} \cong 2\pi/q_{\mathrm{max}} = 2 \, \mathrm{nm}$ and $r_{\mathrm{max}} \cong 2\pi/q_{\mathrm{min}} = 130 \, \mathrm{nm}$.}
\end{figure}

\section{Conclusion}

Using magnetic small-angle neutron scattering (SANS) we have provided a comparative study of the magnetic microstructure of textured and isotropic $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}/\alpha$-$\mathrm{Fe}$ nanocomposites. Our neutron-data analysis suggests that the spin-misalignment scattering of the textured sample is dominated by spin components along one direction perpendicular to the easy $c$-axis (pressing direction) of the $\mathrm{Nd}_2\mathrm{Fe}_{14}\mathrm{B}$ grains. This anisotropy in the magnetization distribution is accompanied by the presence of a crystallographic texture along these directions. Possible origins for the spin canting (on an intraparticle length scale) have been discussed and are related to the presence of perturbed interface regions, crystalline imperfections, and/or a dispersion in the orientation distribution of the easy $c$-axes. In agreement with the x-ray synchrotron and neutron data, we find anisotropic real-space correlations, with a correlation length that has been estimated at about $40-50 \, \mathrm{nm}$ in the remanent state. The results demonstrate the power of magnetic SANS for analyzing anisotropic magnetic structures on a nanometer length scale; in particular, the complimentary use of the perpendicular and parallel scattering geometries has (for the textured sample) provided results that were otherwise not accessible with only one geometry.

\section*{Acknowledgements}

Denis Mettus acknowledges financial support from the National Research Fund of Luxembourg (INTER/DFG/12/07). This paper is based on results obtained from the future pioneering program ``Development of magnetic material technology for high-efficiency motors'' commissioned by the New Energy and Industrial Technology Development Organization (NEDO). The neutron experiments were performed at the Swiss spallation neutron source SINQ, Paul Scherrer Institute, Villigen, Switzerland. ALBA synchrotron is acknowledged for the provision of beamtime. We thank Birgit Heiland (INM, Saarbr\"ucken) and J\"org Schmauch (Universit\"at des Saarlandes) for the electron-microscopy work.

\bibliographystyle{apsrev4-1}