$L^{X^{x}}(\mathrm{d}u,\mathrm{d}z)$ stands for the integration in space and
time with respect to the local time of $X^{x}$ and $v^{\prime}$ is the spatial
partial derivative of $v$ solution to PDE (5.28).

###### Proof.

We recall that, Theorem 3.4 ensures the unique weak solvability of FBSDE
(5.41) and there is a ”decoupling field” $v\in
W^{1,2,d+1}_{\text{loc}}([0,T[\times\mathbb{R}^{d},\mathbb{R})$ solution to
PDE (5.28) such that the forward equation in (5.27) can be written as:

$\displaystyle
X_{t}^{x}=x+\int_{s}^{t}\tilde{b}(r,X_{r}^{x})\mathrm{d}r+W_{t}-W_{s},$ (5.47)

where $\tilde{b}(t,x):=b(t,x,v(t,x),\nabla_{x}v(t,x)).$ Using (3.6) , (3.10)
and the assumptions of Proposition 5.22 the drift $\tilde{b}$ is bounded for
all $t\in[0,T-\delta]$. Thus, using [40, Theorem 3] the forward process
$X^{x}$ is Sobolev differentiable. On the other hand, the function $v$ is
Lipschitz continuous for all $t<T$ (see (3.10) ) then for all $\omega$ outside
every $\mathbb{P}\text{-null}$ set, the process
$Y^{x}_{t}(\omega)=v(t,X_{t}^{x}(\omega))$ belongs to
$W^{1,1}_{Loc}(\mathbb{R})$ (see [33]). Then (5.46) follows by applying the
chain rule developed in the same paper [33, Theorem 1.1]. This concludes the
proof. ∎

######  Remark 5.24.

The following relation between the Sobolev and the Malliavin derivatives of
the forward and backward components hold respectively
$\mathrm{d}t\otimes\mathbb{P}\text{-a.s.}$

$\displaystyle D_{s}X_{t}^{x}$ $\displaystyle=\xi_{t}^{x}(\xi_{s}^{x})^{-1},$
$\displaystyle D_{s}Y_{t}^{x}$ $\displaystyle=\eta_{t}^{x}(\xi_{s}^{x})^{-1}$

for any $s,t\in[0,T),s\leq t.$

####  5.3.2. Density analysis of FBSDEs with discontinuous drift case

The main objective in this section is to conduct the analysis of densities
obtained in the preceding section by assuming now

######  Assumption 5.25.

  * (A9)

Assumption 3.2 is valid in one dimension $(d=1)$ with $\sigma\equiv 1$ and $b$
bounded in the forward and the control variable i.e.,
$\forall(t,x,y,z)\in[0,T]\times\mathbb{R}\times\mathbb{R}\times\mathbb{R},|b(t,x,y,z)|\leq\Lambda(1+|y|)$.

  * (A10)

Assumptions (A4) and (A5) hold.

######  Theorem 5.26 (Existence of a Density for $X_{t}$).

Under (A9), the finite dimensional laws of the strong solution to the forward
SDE (5.47) are absolutely continuous with respect to the Lebesgue measure.
Precisely for fixed $t\in(0,T],$ the law of $X_{t}^{x}$ denoted by
$\mathcal{L}(X_{t}^{x})$ has a density with respect to the Lebesgue measure.

###### Proof.

This follows by showing that the exponential term in equation (5.42) is finite
almost surely. The latter is obtained by using Girsanov theorem and Lemma A.1.
This ends the proof. ∎

#### Existence of a Density for the backward component $Y_{t}$

We first derive the linear BSDE satisfied by a version
$(D_{s}Y_{t}^{x},D_{s}Z_{t}^{x})$ of the Malliavin derivatives of the solution
to (5.41).

######  Proposition 5.27.

Under Assumption 5.25, a version of $(D_{s}Y_{t},D_{s}Z_{t})$ satisfies

$\displaystyle D_{s}Y_{t}$ $\displaystyle=0,\quad D_{s}Y_{t}=0,\quad t<s\leq
T,$ $\displaystyle D_{s}Y_{t}$
$\displaystyle=\phi^{\prime}(X_{T})D_{s}X_{T}+\int_{t}^{T}\langle\nabla
f(r,\mathbb{X}_{r}),D_{s}\mathbb{X}_{r}\rangle\mathrm{d}r-\int_{t}^{T}D_{s}Z_{r}\mathrm{d}W_{r},$
(5.48)

where $D_{s}X_{t}$ is given by (5.42). Moreover, $(D_{t}Y_{t})_{0\leq t\leq
T}$ is a continuous version of $(Z_{t})_{0\leq t\leq T}$.

###### Proof.

The first part of the proof follows from Theorem 5.22. As for the
representation, we appeal [23, Theorem 5.3] since $\tilde{b}$ is uniformly
bounded for all $t\in[0,T]$. The proof is completed. ∎

Below is stated the main result of this subsection.

######  Theorem 5.28.

Suppose Assumption 5.25 is in force. Suppose in addition there is
$\mathcal{A}\in\mathcal{B}(\mathbb{R})$ such that
$\mathbb{P}(X_{T}\in\mathcal{A}|\mathfrak{F}_{t})>0$ and one of the following
assumptions holds

  * (H+)

$\phi^{\prime}\geq 0,$
$\phi^{\prime}_{|\mathcal{A}}>0,\mathcal{L}(X_{T})\text{-a.e.}$ and
$\underline{f}\geq 0,$

  * (H–)

$\phi^{\prime}\leq 0,$
$\phi^{\prime}_{|\mathcal{A}}<0,\mathcal{L}(X_{T})\text{-a.e.}$ and
$\overline{f}\leq 0.$

Then, $Y_{t}$ possesses an absolute continuous law which with respect to the
Lebesgue measure on $\mathbb{R}$.

###### Proof.

From Assumption 5.25, we deduce that

$\displaystyle|f_{z}^{\prime}(s,\mathbb{X}_{s})|\leq\Lambda\big{(}1+\ell(Y_{s})|Z_{s}|\big{)}\in\mathcal{H}_{\text{BMO}}.$
(5.49)

Let $\mathbb{Q}$ be the probability measure defined by

$\displaystyle\frac{\mathrm{d}\mathbb{Q}}{\mathrm{d}\mathbb{P}}\Big{|}_{\mathfrak{F}_{t}}=\mathcal{E}\Big{(}f_{z}^{\prime}(t,\mathbb{X}_{t})*W\Big{)}_{t}=M_{t}.$
(5.50)

It follows from (5.49) that $M_{t}$ is uniformly integrable. Therefore, taking
the conditioning expectation with respect to $\mathbb{Q}$ on both sides of
(5.27) gives for all $0\leq s\leq t<T$

$\displaystyle
D_{s}Y_{t}=\mathbb{E}^{\mathbb{Q}}\Big{(}\phi^{\prime}(X_{T})D_{s}X_{T}+\int_{t}^{T}\left(f_{x}^{\prime}(r,\mathbb{X}_{r})D_{s}X_{r}+f_{y}^{\prime}(s,\mathbb{X}_{r})D_{s}Y_{r}\right)\mathrm{d}r\Big{|}\mathfrak{F}_{t}\Big{)}$
(5.51)

Using Remark 5.24 and the well known linearisation method, $D_{s}Y_{t}$ can be
rewritten:

$\displaystyle D_{s}Y_{t}$
$\displaystyle=\mathbb{E}^{\mathbb{Q}}\Big{[}e^{\int_{t}^{T}f_{y}^{\prime}(u,\mathbb{X}_{u})\mathrm{d}u}\phi^{\prime}(X_{T})\xi_{T}+\int_{t}^{T}e^{\int_{t}^{r}f_{y}^{\prime}(u,\mathbb{X}_{u})\mathrm{d}u}f_{x}^{\prime}(r,\mathbb{X}_{r})\xi_{r}\mathrm{d}r\Big{|}\mathfrak{F}_{t}\Big{]}(\xi_{s})^{-1}$
$\displaystyle=\mathbb{E}\Big{[}\psi_{T}\phi^{\prime}(X_{T})\xi_{T}+\int_{t}^{T}\psi_{r}f_{x}^{\prime}(r,\mathbb{X}_{r})\xi_{r}\mathrm{d}r\Big{|}\mathfrak{F}_{t}\Big{]}(\psi_{t})^{-1}(\xi_{s})^{-1},$

where the last equality is due to the Bayes’ rule with $\psi$ is given by

$\psi_{t}:=\exp\Big{(}\int_{0}^{t}(f_{y}^{\prime}(u,\mathbb{X}_{u})-\frac{1}{2}|f_{z}^{\prime}(u,\mathbb{X}_{u})|^{2})\mathrm{d}u+\int_{0}^{t}f_{z}^{\prime}(u,\mathbb{X}_{u})\mathrm{d}W_{u}\Big{)},$

and the expression of $\xi$ is given in (5.45). Therefore, from the
computations above, the Malliavin covariance $\Gamma_{Y_{t}}$ of $Y_{t}$ is
given by

$\Gamma_{Y_{t}}=\Big{(}\mathbb{E}\Big{[}\psi_{T}\phi^{\prime}(X_{T})\xi_{T}+\int_{t}^{T}\psi_{r}f_{x}^{\prime}(r,\mathbb{X}_{r})\xi_{r}\mathrm{d}r\Big{|}\mathfrak{F}_{t}\Big{]}\Big{)}^{2}\times(\psi_{t}^{-1})^{2}\int_{0}^{t}(\xi_{s}^{-1})^{2}\mathrm{d}s.$

It is enough to prove that
$\mathbb{E}\Big{[}\psi_{T}\phi^{\prime}(X_{T})\xi_{T}+\int_{t}^{T}\psi_{r}f_{x}^{\prime}(r,\mathbb{X}_{r})\xi_{r}\mathrm{d}r\big{|}\mathfrak{F}_{t}\Big{]}\neq
0$ under assumption (H+) to deduce the desired result. Let us remark that, the
product $\psi_{T}\xi_{T}$ can be rewritten as follow:

$\displaystyle\psi_{T}\xi_{T}$
$\displaystyle=\exp\Big{(}-\int_{0}^{T}\int_{\mathbb{R}}\tilde{b}\left(u,z\right)L^{X^{x}}(\mathrm{d}u,\mathrm{d}z)+\int_{0}^{T}f_{y}^{\prime}(u,\mathbb{X}_{u})\mathrm{d}u\Big{)}$
$\displaystyle\quad\times\exp\Big{(}\int_{0}^{T}f_{z}^{\prime}(u,\mathbb{X}_{u})\mathrm{d}W_{u}-\frac{1}{2}\int_{0}^{T}|f_{z}^{\prime}(u,\mathbb{X}_{u})|^{2}\mathrm{d}u\Big{)}=H_{T}\times
M_{T}$

The Bayes’ rule once more yields

$\displaystyle\mathbb{E}\Big{[}\psi_{T}\phi^{\prime}(X_{T})\xi_{T}+\int_{t}^{T}\psi_{r}f_{x}^{\prime}(r,\mathbb{X}_{r})\xi_{r}\mathrm{d}r\Big{|}\mathfrak{F}_{t}\Big{]}=\mathbb{E}^{\mathbb{Q}}\Big{[}\phi^{\prime}(X_{T})H_{T}+\int_{t}^{T}f_{x}^{\prime}(r,\mathbb{X}_{r})H_{r}\mathrm{d}r\Big{|}\mathfrak{F}_{t}\Big{]}M_{t}$

Using the fact that
$\phi^{\prime}+f_{x}^{\prime}\geq\underline{\phi}+\underline{f}$, we deduce:

$\displaystyle\mathbb{E}\Big{[}\psi_{T}\phi^{\prime}(X_{T})\xi_{T}+\int_{t}^{T}\psi_{r}f_{x}^{\prime}(r,\mathbb{X}_{r})\xi_{r}\mathrm{d}r\Big{|}\mathfrak{F}_{t}\Big{]}\geq$
$\displaystyle\mathbb{E}^{\mathbb{Q}}\Big{[}{1}_{\\{X_{T}\in\mathcal{A}\\}}\Big{(}\underline{\phi}H_{T}+\underline{f}\int_{t}^{T}H_{r}\mathrm{d}r\Big{)}\Big{|}\mathfrak{F}_{t}\Big{]}M_{t}$

Using the bound
$|f_{y}^{\prime}(s,\mathbb{X}_{s})|\leq\Lambda(1+|Z_{s}|^{\alpha})$ and Young
inequality $|z|^{\alpha}\leq c(1+|z|^{2})$, for $\alpha\in(0,1)$ we deduce

$\displaystyle\mathbb{E}\Big{[}\psi_{T}\phi^{\prime}(X_{T})\xi_{T}+\int_{t}^{T}\psi_{r}f_{x}^{\prime}(r,\mathbb{X}_{r})\xi_{r}\mathrm{d}r\Big{|}\mathfrak{F}_{t}\Big{]}\geq$
$\displaystyle\mathbb{E}^{\mathbb{Q}}\Big{[}{1}_{\\{X_{T}\in\mathcal{A}\\}}\Big{(}\underline{\phi}^{\mathcal{A}}e^{-\int_{0}^{T}\int_{\mathbb{R}}\tilde{b}\left(u,z\right)L^{X^{x}}(\mathrm{d}u,\mathrm{d}z)+\int_{0}^{T}f_{y}^{\prime}(u,\mathbb{X}_{u})\mathrm{d}u}$
$\displaystyle+\underline{f}\int_{t}^{T}e^{-\int_{0}^{r}\int_{\mathbb{R}}\tilde{b}\left(u,z\right)L^{X^{x}}(\mathrm{d}u,\mathrm{d}z)}e^{\int_{0}^{r}f_{y}^{\prime}(u,\mathbb{X}_{u})\mathrm{d}u}\mathrm{d}r\Big{)}|\mathfrak{F}_{t}\Big{]}M_{t}$
$\displaystyle\geq$
$\displaystyle\mathbb{E}^{\mathbb{Q}}\Big{[}{1}_{\\{X_{T}\in\mathcal{A}\\}}e^{-\Lambda
T}\Big{(}\underline{\phi}^{\mathcal{A}}e^{-\int_{0}^{T}\int_{\mathbb{R}}\tilde{b}\left(u,z\right)L^{X^{x}}(\mathrm{d}u,\mathrm{d}z)}e^{-\Lambda\int_{0}^{T}|Z_{u}|^{2}\mathrm{d}u}$
$\displaystyle+\underline{f}\int_{t}^{T}e^{-\int_{0}^{r}\int_{\mathbb{R}}\tilde{b}\left(u,z\right)L^{X^{x}}(\mathrm{d}u,\mathrm{d}z)}e^{-\Lambda\int_{0}^{T}|Z_{u}|^{2}\mathrm{d}u}\Big{)}|\mathfrak{F}_{t}\Big{]}M_{t}.$

Provided that the law of $X_{T}$ is absolutely continuous with respect to the
Lebesgue measure, and that
$-\int_{0}^{T}\int_{\mathbb{R}}\tilde{b}\left(u,z\right)L^{X^{x}}(\mathrm{d}u,\mathrm{d}z)<\infty,\,\,\mathbb{Q}$-a.s.,
we deduce that
$\mathbb{E}\Big{(}\psi_{T}\phi^{\prime}(X_{T})\xi_{T}+\int_{t}^{T}\psi_{r}f_{x}^{\prime}(r,\mathbb{X}_{r})\xi_{r}\mathrm{d}r\Big{|}\mathfrak{F}_{t}\Big{)}>0$.
The absolute continuity of the law of with respect to the Lebesgue measure
follows from assumption. It then remains to show that:

$-\int_{0}^{T}\int_{\mathbb{R}}\tilde{b}\left(u,z\right)L^{X^{x}}(\mathrm{d}u,\mathrm{d}z)<\infty,\,\,\mathbb{Q}\text{-a.s.}$

Thanks to Girsanov’s theorem, the decomposition (A) and Hölder inequality we
have

$\displaystyle\mathbb{E}^{\mathbb{Q}}\Big{[}-\int_{0}^{T}\int_{\mathbb{R}}\tilde{b}\left(s,z\right)L^{X^{x}}(\mathrm{d}s,\mathrm{d}z)\Big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}-M_{T}\int_{0}^{T}\int_{\mathbb{R}}\tilde{b}\left(s,z\right)L^{X^{x}}(\mathrm{d}s,\mathrm{d}z)\Big{]}\leq\mathbb{E}\Big{[}M_{T}^{2}\Big{]}^{\frac{1}{2}}\mathbb{E}\Big{[}\Big{(}-\int_{0}^{T}\int_{\mathbb{R}}\tilde{b}\left(s,z\right)L^{X^{x}}(\mathrm{d}s,\mathrm{d}z)\Big{)}^{2}\Big{]}^{\frac{1}{2}}$
$\displaystyle\leq$ $\displaystyle
C\mathbb{E}\Big{[}e^{2\int_{0}^{T}\tilde{b}(s,W_{s}^{x})\mathrm{d}W_{s}-\int_{0}^{T}(\tilde{b}(s,W_{s}^{x}))^{2}\mathrm{d}s}\Big{(}-\int_{0}^{T}\int_{\mathbb{R}}\tilde{b}\left(s,z\right)L^{W^{x}}(\mathrm{d}s,\mathrm{d}z)\Big{)}^{2}\Big{]}^{\frac{1}{2}}$
$\displaystyle\leq$ $\displaystyle
C\mathbb{E}\Big{[}e^{2\int_{0}^{T}\tilde{b}(s,W_{s}^{x})\mathrm{d}W_{s}-\int_{0}^{T}(\tilde{b}(s,W_{s}^{x}))^{2}\mathrm{d}s}\Big{(}-\int_{0}^{T}\tilde{b}(s,W^{x}_{s})\mathrm{d}W_{s}$
$\displaystyle-\int_{0}^{T}\tilde{b}(T-s,\widehat{W}^{x}_{s})\mathrm{d}B(s)+\int_{0}^{T}\tilde{b}(T-s,\widehat{W}^{x}_{s})\frac{\widehat{W}_{s}}{T-s}\mathrm{d}s\Big{)}^{2}\Big{]}^{\frac{1}{2}}$
$\displaystyle\leq$ $\displaystyle
C\mathbb{E}\Big{[}e^{4\int_{0}^{t}\tilde{b}(s,W_{s}^{x})\mathrm{d}W_{s}-2\int_{0}^{t}(\tilde{b}(s,W_{s}^{x}))^{2}\mathrm{d}s}\Big{]}^{1/4}\Big{\\{}\mathbb{E}\Big{[}\Big{(}\int_{0}^{T}\tilde{b}(s,W^{x}_{s})\mathrm{d}W_{s}\Big{)}^{4}\Big{]}^{1/4}$
$\displaystyle+\mathbb{E}\Big{[}\Big{(}\int_{0}^{T}\tilde{b}(T-s,\widehat{W}^{x}_{s})\mathrm{d}B_{s}\Big{)}^{4}\Big{]}^{1/4}+\mathbb{E}\Big{[}\Big{(}\int_{0}^{T}\tilde{b}(T-s,\widehat{W}^{x}_{s})\frac{\widehat{W}_{s}}{T-s}\mathrm{d}s\Big{)}^{4}\Big{]}^{1/4}\Big{\\}}$
$\displaystyle\leq$ $\displaystyle CI_{1}(I_{2}+I_{3}+I_{4}).$ (5.52)

Using Cauchy-Schwartz inequality, we have that

$I_{1}\leq\mathbb{E}\Big{[}e^{4\int_{0}^{t}\tilde{b}(s,W_{s}^{x})\mathrm{d}W_{s}-8\int_{0}^{t}(\tilde{b}(s,W_{s}^{x}))^{2}\mathrm{d}s}\Big{]}^{1/8}\mathbb{E}\Big{[}e^{6\int_{0}^{t}(\tilde{b}(s,W_{s}^{x}))^{2}\mathrm{d}s}\Big{]}^{1/8}.$

Applying Girsanov theorem to the martingale
$4\int_{0}^{t}\tilde{b}(s,W_{s}^{x})\mathrm{d}W_{s}$, the first term in the
bound above is equal to $1$. The boundedness of $\tilde{b}$ yields the
boundedness of the second term.

Using Itô isometry (or Burkholder-Davis-Gundy), the boundedness of $\tilde{b}$
implies the boundedness of $I_{2}$, $I_{3}$ and $I_{4}$ (see for e.g. [7]).
The proof is completed. ∎

######  Remark 5.29.

Another interpretation of the result in Theorem 5.28 is that the spatial
derivative $v^{\prime}$ of $v$ (solution to quasi-linear PDE (5.28)) is sucht
that $v^{\prime}(t,x)>0\text{ a.e.}$ (or $v^{\prime}(t,x)<0\text{ a.e.}$).
This follows from the fact that $D_{t}Y_{t}$ is a continuous version of
$Z_{t}=v^{\prime}(t,x)$ for all $t\in(0,T)$ (see [2, Remark 3.4]). We
emphasize that the solution $v$ here is understood in the weak sense. In other
words, we obtain conditions that guarantee the strict monotonicity of weak
solutions of parabolic PDE by means of purely probabilistic tools. To our
knowledge, such investigation has always been performed only for classical
solution of parabolic or elliptic PDEs.

##  6\. Application to smoothness of densities

The main aim in this section is to apply the results obtained above to study
the regularity of densities of solution to coupled FBSDEs with quadratic
drivers under minimal conditions. Indeed, it is well known that the Malliavin
calculus stands as a powerful tool to show hypoellipticity (existence of
smooth density) under weaker conditions than Hörmander’s ones (see [41]). In
particular, we exploit Theorem 5.13 to prove that the densities obtained for
$X$, $Y$ and $Z$ in the above sections, are indeed Hölder continuous despite
the fact that the drift coefficient is at most only Lipschitz continuous.

We consider the following FBSDE

$\begin{cases}X_{t}=x+\displaystyle\int_{0}^{t}b(s,X_{s},Y_{s},Z_{s})\mathrm{d}s+W_{t},\\\
Y_{t}=\phi(X_{T})+\displaystyle\int_{t}^{T}f(s,X_{s},Y_{s},Z_{s})\mathrm{d}s-\int_{t}^{T}Z_{s}\mathrm{d}W_{s}.\end{cases}$
(6.1)

Recall that under assumption (A6) and if $\phi^{\prime}\geq 0$,
$\underline{f}\geq 0$ (resp. $\phi^{\prime}\leq 0$ and $\underline{f}\leq 0$)
then for all $0<s\leq t<T,$ $D_{s}Y_{t}\geq 0$ (resp. $D_{s}Y_{t}\leq 0$ ).
Moreover, if there is a Borel set $\mathcal{A}$ such that
$\mathbb{P}(X_{T}\in\mathcal{A}/\mathfrak{F}_{t})>0$ and
$\phi^{\prime}_{|\mathcal{A}}>0$ (resp. $\phi^{\prime}_{|\mathcal{A}}<0$) then
$D_{s}Y_{t}>0$ (resp. $D_{s}Y_{t}<0$ ). In addition, if (A7) and (A8) are also
valid, then the processes $X,Y$ and $Z$ solutions to (6.1) are twice Malliavin
differentiable.

Let $(G_{k})_{0\leq k\leq n}$ be a random variable such that $G_{0}=1$ and for
all $0\leq s\leq t\leq T$

$G_{k}=\delta\Big{(}G_{k-1}\Big{(}\int_{0}^{T}D_{s}X_{t}\mathrm{d}s\Big{)}^{-1}\Big{)},\text{
if }1\leq k\leq n+1,$

with $\delta$ denoting the Skorohod integral introduced in (2.1). From [41,
P.115], if $X_{t}\in\mathbb{D}^{1,2}(\Omega)$,
$\int_{0}^{T}D_{s}X_{t}\mathrm{d}s\neq 0$ $\mathbb{P}\text{-a.s.}$ and
$G_{k}\Big{(}\int_{0}^{T}D_{s}X_{t}\mathrm{d}s\Big{)}^{-1}$ belongs to
$\text{Dom}(\delta)$ then the probability density $\rho_{X_{t}}$ of the
forward equation $X_{t}$ is of class $C^{n}(\mathbb{R})$ and

$\displaystyle\frac{d^{n}}{dx^{n}}\rho_{X_{t}}(x)=(-1)^{n}\mathbb{E}\left[1_{\\{X_{t}>x\\}}G_{n+1}\right]$
(6.2)

In the above context, the Malliavin derivative of $X_{t}$ is explicitly given
by

$D_{s}X_{t}=\exp\Big{(}\int_{s}^{t}\tilde{b}^{\prime}(u,X_{u})\mathrm{d}u\Big{)},\text{
for all }0\leq s\leq t,$

where $\tilde{b}(t,x)=b(t,x,v(t,x),v^{\prime}(t,x))$ and $\tilde{b}^{\prime}$
denotes the weak derivative of $\tilde{b}$ with respect to the variable $x$.
Hence for any $t\in(0,T]$ there is $\epsilon>0$ such that
$\int_{0}^{T}D_{s}X_{t}\mathrm{d}s\geq\epsilon>0.$ Using the smoothness of the
function $y\mapsto y^{-1}$ on $(\epsilon,\infty)$ and the fact that
$X_{t}\in\mathbb{D}^{2,2}(\Omega)$ we deduce that
$F:=\Big{(}\int_{0}^{T}D_{s}X_{t}\mathrm{d}s\Big{)}^{-1}\in\mathbb{D}^{1,2}(\Omega)$,
so $F\in\text{Dom}(\delta)$ and $G_{1}=FW(T)-\int_{0}^{T}D_{s}F\mathrm{d}s$
which implies that $G_{1}\in\mathbb{D}^{1,2}(\Omega)$, then
$G_{1}F\in\mathbb{D}^{1,2}(\Omega)$ and $G_{1}F\in\text{Dom}(\delta)$.

Similarly, let $(H_{k})_{0\leq k\leq n}$ be a random variable such that
$H_{0}=1$ and for all $0\leq s\leq t\leq T$

$H_{k}=\delta\Big{(}H_{k-1}\Big{(}\int_{0}^{T}D_{s}Y_{t}\mathrm{d}s\Big{)}^{-1}\Big{)},\text{
if }1\leq k\leq n+1.$

Since for all $0\leq s\leq t$, $D_{s}Y_{t}>0$ $\mathbb{P}\text{-a.s.}$, we
deduce that, there is a constant $\epsilon_{1}>0$ such that
$\int_{0}^{T}D_{s}Y_{t}\mathrm{d}s\geq\epsilon_{1}$. Using now the smoothness
of the function $y\mapsto y^{-1}$ on $(\epsilon_{1},\infty)$ and the fact that
$Y_{t}\in\mathbb{D}^{2,2}(\Omega)$ we deduce that
$F_{1}:=\Big{(}\int_{0}^{T}D_{s}Y_{t}\mathrm{d}s\Big{)}^{-1}\in\mathbb{D}^{1,2}(\Omega)$
so $H_{1}F_{1}\in\text{Dom}(\delta)$.

We have shown that

######  Proposition 6.1.

Under (A6),(A7),(A8) and assuming further that either (A+) or (A-) is valid.
Then, the densities $\rho_{X_{t}}$ and $\rho_{Y_{t}}$ of the processes $X_{t}$
and $Y_{t}$ solution to (6.1), respectively belong to $C^{0}(\mathbb{R})$ such
that

$\displaystyle\rho_{X_{t}}(x)$
$\displaystyle=\mathbb{E}\left(1_{\\{X_{t}>x\\}}G_{1}\right)$
$\displaystyle\rho_{Y_{t}}(x)$
$\displaystyle=\mathbb{E}\left(1_{\\{Y_{t}>x\\}}H_{1}\right).$

Moreover, if the assumptions of Theorem 5.19 are valid then for all $0\leq
s\leq t$, $D_{s}Z_{t}>0$ $\mathbb{P}\text{-a.s.}$ and the density
$\rho_{Z_{t}}$ of the solution process $Z$ to (6.1) belongs to
$C^{0}(\mathbb{R})$ and

$\displaystyle\rho_{Z_{t}}(x)=\mathbb{E}\left(1_{\\{Z_{t}>x\\}}I_{1}\right),$

where
$I_{1}=\delta\Big{(}\Big{(}\int_{0}^{T}D_{s}Z_{t}\mathrm{d}s\Big{)}^{-1}\Big{)}$.

The following result improves the regularity of the densities of solutions
$X,Y$ and $Z$ to the system (6.1) obtained above. The method consists to
investigate the integrability property of the Malliavin covariance matrix of
the different processes (see Lemma 6.2 below).

We recall that for any random variable $F$, the matrix
$\Gamma_{F}=(\Gamma_{F}^{ij})_{i,j=1,\ldots,d}$, given by
$\Gamma_{F}^{ij}:=\langle DF^{i},DF^{j}\rangle_{\mathcal{H}^{2}}$ is said to
be non-degenerate if for all $p\geq 1$

$\mathbb{E}\Big{[}\left(\Gamma_{F}\right)^{-p}\Big{]}<\infty.$ (6.3)

######  Lemma 6.2 (Proposition 23 in [5]).

Let $F=(F^{1},\cdots,F^{d})$ with $F^{1},\cdots,F^{d}\in\bigcap_{p\geq
1}\mathbb{D}^{k+1,p}$, such that $\Gamma_{F}$ satisfies (6.3). Then
$\rho_{F}\in C^{k-1,\beta}(\mathbb{R})$ for some $k\geq 1$ and $\beta\in(0,1)$
i.e. $\rho_{F}$ is $k-1\text{-times}$ differentiable with Hölder continuous
derivatives of exponent $\beta<1$.

We will first prove that the densities of $X_{t}$ and $Y_{t}$ are Hölder
continuous.

######  Proposition 6.3.

Let assumption (A6),(A7) and (A8) be in force. Suppose in addition that there
is $\mathcal{A}\in\mathcal{B}(\mathbb{R})$ such that
$\mathbb{P}(X_{T}\in\mathcal{A}|\mathfrak{F}_{t})>0$ and one of the
assumptions (A+) or (A–) holds. Then, $\rho_{X_{t}}\in
C^{0,\beta}(\mathbb{R})$ for all $t\in[0,T]$ and for all $t\in[0,T)$, $Y_{t}$
has a $\beta\text{-H\"{o}lder}$ continuous density $\rho_{Y_{t}}$ given by

$\rho_{Y_{t}}(x)=\rho_{Y_{t}}(x_{0})\exp\Big{(}-\int_{x_{0}}^{x}w_{Y_{t}}(z)\mathrm{d}z\Big{)},\,x\in\operatorname{supp}\rho_{Y_{t}},$

where $x_{0}$ is a point in the interior of
$\operatorname{supp}(\rho_{Y_{t}})$ and

$w_{Y_{t}}(z):=\mathbb{E}\Big{[}\delta\Big{(}\int_{0}^{T}D_{s}Y_{t}\mathrm{d}s\Big{)}^{-1}\Big{|}Y_{t}=z\Big{]},$

$\delta$ denotes the Skorohod integral introduced in (2.1).

###### Proof.

From Lemma 5.9, we know that $\tilde{b}$ is uniformly bounded and Lipschitz
continuous in its spatial variable for all $t\in[0,T)$. Then, using [8,
Proposition 4.4] the following holds

$\displaystyle\Big{(}\det\Gamma_{X_{t}}\Big{)}^{-1}\in\bigcap_{p\geq
1}L^{p}(\Omega).$

This implies that $\rho_{X_{t}}\in C^{0,\beta}(\mathbb{R})$ for all
$t\in[0,T)$ with $\beta\in(0,1)$. On the other hand, for all $t\in[0,T)$ we
deduce

$\displaystyle\det\Gamma_{Y_{t}}=(v^{\prime}(t,X_{t}))^{2}\det\Gamma_{X_{t}}\geq\alpha^{2}(t)\det\Gamma_{X_{t}},$

where $\alpha(t)$ is the positive function such that
$v^{\prime}(t,X_{t})\geq\alpha(t)$ or $v^{\prime}(t,X_{t})\leq-\alpha(t)$
$\mathbb{P}\text{-a.s.}$ (see proof of Theorem 5.11 ). Without loss of
generality, we assume that there is a constant $\epsilon_{1}>0$ such that
$\alpha(t)\geq\epsilon_{1}$. Hence,
$\Big{(}\det\Gamma_{Y_{t}}\Big{)}^{-1}\leq\alpha^{-2}(t)\Big{(}\det\Gamma_{X_{t}}\Big{)}^{-1}\in\bigcap_{p\geq
1}L^{p}(\Omega)$. For the explicit representation of the density it follows by
combining arguments of the proof of [41, Proposition 2.1.1] and [15,
Proposition 2]. This concludes the proof. ∎

Let us turn now to the case of the control process $Z$.

######  Proposition 6.4.

Assume that there is $\mathcal{A}\in\mathcal{B}(\mathbb{R})$ such that
$\mathbb{P}(X_{T}\in\mathcal{A}/\mathfrak{F}_{t})>0$ and let assumptions of
Theorem 5.19 be in force. Then the control process $Z$ solution to (6.1) has a
density $\rho_{Z_{t}}$ which is Hölder continuous and the following explicit
representation holds

$\rho_{Z_{t}}(z)=\rho_{Z_{t}}(x_{0})\exp\Big{(}-\int_{z_{0}}^{z}w_{Z_{t}}(u)\mathrm{d}u\Big{)},\,z\in\operatorname{supp}\rho_{Z_{t}},$

where $z_{0}$ is a point in the interior of
$\operatorname{supp}(\rho_{Z_{t}})$ and

$w_{Z_{t}}(u):=\mathbb{E}\Big{[}\delta\Big{(}\int_{0}^{T}D_{s}Z_{t}\mathrm{d}s\Big{)}^{-1}\Big{|}Z_{t}=u\Big{]},$

$\delta$ denotes the Skorohod integral introduced in (2.1).

###### Proof.

We recall that for all $t\in[0,T)$ $Z_{t}=v^{\prime}(t,X_{t})$
$\mathbb{P}\text{-a.s.}$ Hence

$\displaystyle\det\Gamma_{Z_{t}}=(v^{\prime\prime}(t,X_{t}))^{2}\det\Gamma_{X_{t}}\geq\omega(t)^{2}\det\Gamma_{X_{t}},$

where $\omega(t)$ is the positive function such that
$v^{\prime\prime}(t,X_{t})\geq\omega(t)$ or
$v^{\prime\prime}(t,X_{t})\leq-\omega(t)$ $\mathbb{P}\text{-a.s.}$ (see proof
of Corollary 5.20 ). Then, for all $p\geq 1$

$\mathbb{E}\left[\Big{(}\det\Gamma_{Z_{t}}\Big{)}^{-p}\right]\leq\omega^{-2p}(t)\mathbb{E}\left[\Big{(}\det\Gamma_{X_{t}}\Big{)}^{-p}\right]<\infty.$

∎

### Conclusion and Possible extension

In this paper, we derive some conditions under which the solution of a system
of coupled forward-backward SDEs with non-smooth drift coefficient and
quadratic driver admits a density which is absolutely continuous with respect
to the Lebesgue measure. Two cases were considered: the Hölder continuous
drift case and the discontinuous drift one. In the latter, we only consider a
Brownian motion with a drift coefficient which in fact depends on the solution
$(Y,Z)$ of the backward SDE. It is worth noting that, there is less hope to
obtaining such densities’ results on coupled FBSDEs with non constant
diffusive coefficient while the drift coefficient is allowed to be
discontinuous. Indeed, via the decoupling field the forward equation reduces
to a classical SDE with singular drifts and the solvability of such equations
with non constant volatility remains a big challenge.

Due to the quadratic behavior of the driver, most of the results derived in
this paper are restricted to the one dimension setting. However, the result on
the weak solvability of FBSDEs obtained in Section 3 (see Theorem 3.4) can be
extended in the multidimensional setting if one considers that the driver $f$
is at most linear growth in $z$ (see Section 8 in [12]). In this case, the
solution $(X,Y,Z)$ will be taken in the space
$\mathcal{S}^{2}(\mathbb{R}^{d})\times\mathcal{S}^{2}(\mathbb{R}^{q})\times\mathcal{H}^{2}(\mathbb{R}^{q\times
d})$, where $q$ represents the dimension of the backward component $Y$. In
this context, the strong solvability results of FBSDEs with the additional
assumptions that the coefficients are Hölder continuous derived in Section 5
remain valid as well as the density’s result of the forward component $X$. The
case for the backward component $Y$ requires more care, since the question
relative to the strict positivity of solution to multidimensional quasi-linear
PDE remains a non trivial problem to handle. However, by combining the
techniques developed in this paper and the recent significant development in
[11] we believe that one can successively handle this question. Many other
interesting generalisation of this paper are conceivable, among them we point
out a possible extension to the multidimensional quadratic setting of fully
coupled FBSDEs by using the recent results from [9] and [25]. These questions
are under investigation for future work.

##  Appendix A Integration with local time

Here, we provide some basic notions related to the integration with local
time. We refer the reader to [17], [18] and [7] for more information. The
integration with respect to local time-space is defined for integrands that
taking value in the space $(\mathcal{H}_{x},\|\cdot\|_{x})$ (see e.g. [17])
representing the space of Borel measurable functions
$\phi:[0,T]\times\mathbb{R}\rightarrow\mathbb{R}$ with the norm

$\displaystyle\left\|\phi\right\|_{x}$
$\displaystyle:=2\Big{(}\int_{0}^{1}\int_{\mathbb{R}}\phi^{2}(s,z)\exp(-\frac{|z-x|^{2}}{2s})\frac{\mathrm{d}s\,\mathrm{d}z}{\sqrt{2\pi
s}}\Big{)}^{\frac{1}{2}}+\int_{0}^{1}\int_{\mathbb{R}}|z-x||\phi(s,x)|\exp(-\frac{|z-x|^{2}}{2s})\frac{\mathrm{d}s\,\mathrm{d}z}{s\sqrt{2\pi
s}}.$

In addition, any bounded function $\phi$ belongs to $\mathcal{H}_{x}$ and the
following representation holds for all $t\in[0,T]$ (see for instance [7, Lemma
2.11]) :

$\displaystyle\int_{0}^{t}\partial_{x}\phi(s,X^{x}(s))\mathrm{d}s=-\int_{0}^{t}\int_{\mathbb{R}}\phi(s,z)L^{X^{x}}(\mathrm{d}s,\mathrm{d}z).$
(A.1)

Furthermore for $\phi\in{\mathcal{H}}_{0}$, the subsequent decomposition is
valid (see for example [18, Theorem 2.1])

$\displaystyle\int_{0}^{t}\int_{\mathbb{R}}\phi(s,z)L^{W^{x}}(\mathrm{d}s,\mathrm{d}z)$
$\displaystyle=\int_{0}^{t}\phi(s,W^{x}_{s})\mathrm{d}W_{s}+\int_{T-t}^{T}\phi(T-s,\widehat{W}^{x}(s))\mathrm{d}B(s)$
$\displaystyle\quad-\int_{T-t}^{T}\phi(T-s,\widehat{W}^{x}(s))\frac{\widehat{W}(s)}{T-s}\mathrm{d}s,\quad
0\leq t\leq T\text{ a.s.},$ (A.2)

where $L^{W^{x}}(\mathrm{d}s,\mathrm{d}z)$ denotes integration with respect to
the local time of the Brownian motion $W^{x}$ in both time and space,
$W^{x}_{\cdot}:=x+W_{\cdot}$ is the Brownian motion started at $x$ and
$\widehat{W}$ is the time-reversed Brownian motion defined by

$\displaystyle\widehat{W}_{t}:=W_{T-t},\,\,0\leq t\leq T.$ (A.3)

The process $B=\\{B(t),\,\,\,0\leq t\leq T\\}$ stands for an independent
Brownian motion with respect to the filtration $\mathcal{F}_{t}^{\widehat{W}}$
generated by $\widehat{W}_{t}$, and satisfies:

$\displaystyle
B_{t}=\widehat{W}_{t}-W_{T}+\int_{t}^{T}\frac{\widehat{W}_{s}}{T-s}\mathrm{d}s.$
(A.4)

We borrow the following result in [7, Lemma A.2]. It gives us the exponential
bound of the local time-space integral of any bounded function.

######  Lemma A.1.

Let $b:[0,T]\times\mathbb{R}\rightarrow\mathbb{R}$ be a bounded and measurable
function. Then for $t\in[0,T],\,\lambda\in\mathbb{R}$ and compact subset
$K\subset\mathbb{R}$, we have

$\underset{x\in
K}{\sup}E\Big{[}\exp\Big{(}\lambda\int_{0}^{t}\int_{\mathbb{R}}b(s,y)L^{W^{x}}(\mathrm{d}s,\mathrm{d}y)\Big{)}\Big{]}<C(\|b\|_{\infty}),$

where $C$ is an increasing function and $L^{W^{x}}(\mathrm{d}s,\mathrm{d}y)$
denotes integration with respect to the local time of the Brownian motion
$W^{x}$ in both time and space.

##  Appendix B Comparison theorem for quadratic BSDE

Consider a BSDE of the form

$Y_{t}=\xi+\int_{t}^{T}f(s,Y_{s},Z_{s})\mathrm{d}s-\int_{t}^{T}Z_{s}\mathrm{d}W_{s},$
(B.1)

satisfying the following set of assumptions

######  Assumption B.1.

  * There exist nonegative constant $\Lambda,K$ and a positive locally bounded function $\ell:\mathbb{R}\mapsto\mathbb{R}_{+}$ such that $\mathbb{P}\text{-a.s.}$

  * (i)

$\xi$ is an $\mathfrak{F}_{T}$ measurable uniformly bounded random variable
i.e.,$\|\xi\|_{L^{\infty}}\leq\Lambda.$

  * (ii)

$f$ is an $\mathfrak{F}$ predictable continuous function and for all
$(\omega,t,y,z)\in\Omega\times[0,T]\times\mathbb{R}\times\mathbb{R}^{d},$

$\displaystyle|f(\omega,t,y,z)|\leq\Lambda(1+|y|+\ell(y)|z|^{2}).$

  * (iii)

For all
$(\omega,t,y,z),(\omega,t,y^{\prime},z^{\prime})\in\Omega\times[0,T]\times\mathbb{R}\times\mathbb{R}^{d},$

$\displaystyle|f(\omega,t,y,z)-f(\omega,t,y^{\prime},z^{\prime})|$
$\displaystyle\leq
K\Big{(}1+\ell(|y-y^{\prime}|^{2})(|z|+|z^{\prime}|)\Big{)}\Big{(}|y-y^{\prime}|+|z-z^{\prime}|\Big{)}.$

Here is the main Theorem developed in this section. We mimic the proof of
Theorem 3.12 in [22].

######  Theorem B.2.

Let $(Y^{i},Z^{i})\in\mathcal{S}^{\infty}\times\mathcal{H}_{\text{BMO}}$ be
the solution to BSDE(B.1), with terminal value $\xi^{i}$ and generator $g^{i}$
for $i\in\\{1,2\\}.$ Under Assumption(B.1) and that

$\displaystyle\xi^{1}\leq\xi^{2},\text{ and }f^{1}(t,Y_{t}^{2},Z_{t}^{2})\leq
f^{2}(t,Y_{t}^{2},Z_{t}^{2})\quad\mathrm{d}\text{t}\otimes\mathrm{d}\mathbb{P}\text{-a.s.}.$

Then for all $t\in[0,T]$ $Y^{1}_{t}\leq Y^{2}_{t}\quad\mathbb{P}\text{-a.s.}.$

Moreover, if either $\xi^{1}<\xi^{2}$ or
$f^{1}(t,Y_{t}^{2},Z_{t}^{2})<f^{2}(t,Y_{t}^{2},Z_{t}^{2})$ in a set of
positive $\mathrm{d}\text{t}\otimes\mathrm{d}\mathbb{P}\text{-measure}$ then,
$Y_{0}^{1}<Y_{0}^{2}.$

###### Proof.

Let us first recall that $\mathcal{S}^{\infty}\text{-norms of }Y^{1},Y^{2}$
are bounded by $\Lambda$ (see e.g. [3, Theorem 3.1]), and the
$\mathcal{H}_{\text{BMO}}\text{-norms of }Z^{1},Z^{2}$ are bounded by a
constant $C$ only depending on $\Lambda,$ and $T$ ( see [3, Corollary 3.2 ]).

Se

$\displaystyle\begin{cases}\delta Y=Y^{1}-Y^{2},&\delta Z=Z^{1}-Z^{2}\\\
\delta\xi=\xi^{1}-\xi^{2},&\delta
f=f^{1}(\cdot,Y^{2},Z^{2})-f^{2}(\cdot,Y^{2},Z^{2}),\end{cases}$

and define the processes

$\displaystyle\Gamma_{t}:$
$\displaystyle=\frac{f^{1}(t,Y_{t}^{1},Z_{t}^{1})-f^{1}(t,Y_{t}^{2},Z_{t}^{1})}{Y^{1}_{t}-Y^{2}_{t}}{1}_{\\{Y^{1}_{t}-Y^{2}_{t}\neq
0\\}},\quad e_{t}:=\exp\Big{(}\int_{0}^{t}\Gamma_{s}\mathrm{d}s\Big{)}$ (B.2)
$\displaystyle\Pi_{t}:$
$\displaystyle=\frac{f^{1}(t,Y_{t}^{2},Z_{t}^{1})-f^{1}(t,Y_{t}^{2},Z_{t}^{2})}{|Z_{t}^{1}-Z_{t}^{2}|^{2}}(Z^{1}_{t}-Z^{2}_{t}){1}_{\\{|Z^{1}_{t}-Z^{2}_{t}|\neq
0\\}}.$ (B.3)

Observe that $|\Pi|\leq\Lambda\left(1+\ell(0)(|Z^{1}|+|Z^{2}|)\right)$, thus
$\int_{0}^{\cdot}|\Pi_{s}|\mathrm{d}B_{s}$ is a BMO martingale, since
$Z^{1},Z^{2}\in\mathcal{H}_{\text{BMO}}.$ Therefore, the probability measure
$\tilde{\mathbb{P}}$ given by
$\mathrm{d}\tilde{\mathbb{P}}/\mathrm{d}\mathbb{P}=\mathcal{E}(\int_{0}^{\cdot}\Pi_{s}\mathrm{d}W_{s})$
is well defined (see [28]) and the process defined by
$W_{\cdot}^{\tilde{\mathbb{P}}}=W_{\cdot}-\int_{0}^{\cdot}\Pi_{s}\mathrm{d}s$
is a $\tilde{\mathbb{P}}\text{-Brownian motion}.$ We deduce from [28], the
existence of $p^{*}>1$ such that
$\mathcal{E}(\int_{0}^{\cdot}\Pi_{s}\mathrm{d}W_{s})\in L^{p^{*}}.$

Since $Y^{1},Y^{2}$ are bounded and $\ell$ is locally bounded, we have
$|\Gamma_{t}|\leq\Lambda(1+2\ell(|Y_{t}^{1}-Y_{t}^{2}|^{2})|Z_{t}^{1}|)\leq\Lambda(1+|Z_{t}^{1}|)$.
Using the same arguments as before, we have
$\Gamma\in\mathcal{H}_{\text{BMO}},$ thus $\Gamma\in\mathcal{H}^{2p}$ for
every $p\geq 1.$

The dynamics of $(\delta Y_{t})$ is given by:

$\displaystyle\delta Y_{t}=\delta\xi+\int_{t}^{T}\delta
f_{s}\mathrm{d}s+\int_{t}^{T}\left(f^{1}(s,Y_{s}^{1},Z_{s}^{1})-f^{1}(s,Y_{s}^{2},Z_{s}^{2})\right)\mathrm{d}s-\int_{t}^{T}\delta
Z_{s}\mathrm{d}W_{s}.$ (B.4)

The Itô’s formula yields

$\displaystyle e_{t}\delta Y_{t}=e_{t}\delta\xi+\int_{t}^{T}e_{s}\delta
f_{s}\mathrm{d}s-\int_{t}^{T}e_{s}\delta
Z_{s}\mathrm{d}W_{s}^{\tilde{\mathbb{P}}},$ (B.5)

where $\int_{0}^{\cdot}e_{s}\delta Z_{s}\mathrm{d}W_{s}^{\tilde{\mathbb{P}}}$
is a true $\tilde{\mathbb{P}}\text{-martingale}.$ Indeed,

$\displaystyle\mathbb{E}^{\tilde{\mathbb{P}}}\int_{0}^{T}|e_{s}|^{2}|\delta
Z_{s}|^{2}\mathrm{d}s$
$\displaystyle\leq\mathbb{E}\Big{[}\Big{(}\mathcal{E}(\int_{0}^{\cdot}\Pi_{s}\mathrm{d}W_{s})\Big{)}^{p^{*}}\Big{]}^{\frac{1}{p^{*}}}\mathbb{E}\Big{[}|e_{T}|^{2q^{*}}\Big{(}\int_{0}^{T}|\delta
Z_{s}|^{2}\mathrm{d}s\Big{)}{q^{*}}\Big{]}^{\frac{1}{p^{*}}},$

and this is finite thanks to the $\mathcal{H}_{\text{BMO}}$ property of
$\Pi,\Gamma,Z^{1}$ and $Z^{2}$. Then, equation (B.5) becomes

$\displaystyle e_{t}\delta
Y_{t}=\mathbb{E}^{\tilde{\mathbb{P}}}\Big{[}e_{t}\delta\xi+\int_{t}^{T}e_{s}\delta
g_{s}\mathrm{d}s|\mathfrak{F}_{t}\Big{]}.$ (B.6)

Thanks to the Theorem’s assumption we conclude that $e_{t}\delta Y_{t}\leq 0$
and hence for all $t\in[0,T]$ we have $Y^{1}_{t}\leq
Y^{2}_{t}\quad\tilde{\mathbb{P}}\text{-a.s.}$ and $\mathbb{P}\text{-a.s.}.$

In particular, for $t=0$ and if $\delta<0$ or if $\delta g<0$ in a set of
positive $\mathrm{d}\text{t}\otimes\mathrm{d}\mathbb{P}\text{-measure},$ then
we conclude that $Y_{0}^{1}<Y_{0}^{2},$ and this ends the proof. ∎

##  Appendix C Second order Malliavin differentiability

In this appendix, we provide some conditions under which a BSDE admits a
second order Malliavin differentiable solution. This result is taken from
[24]. Let us consider the following BSDEs (which can be seen as the BSDE
satisfied by the Malliavin derivative of the backward process )

$\displaystyle U_{s,t}$ $\displaystyle=0,\quad V_{s,t}=0,\quad t\in[0,s),$
$\displaystyle U_{s,t}$
$\displaystyle=D_{s}\xi+\int_{t}^{T}G(r,s,U_{s,r},V_{s,r})\mathrm{d}r-\int_{t}^{T}V_{s,r}\mathrm{d}W_{r},\quad
t\in[s,T]$ (C.1)

where $h$ is a measurable function
$\Omega\times[0,T]\times[0,T]\times\mathbb{R}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$
given by

$G(\cdot,t,s,y,z)=A_{s,t}(\cdot)+B_{t}(\cdot)y+\langle C_{t}(\cdot),z\rangle$

and $\xi$ is an $\mathfrak{F}_{T}\text{-measurable}$ random variable,
satisfying the following assumptions:

  * (H1)

$\xi$ is Malliavin differentiable, $\mathfrak{F}_{T}\text{-measurable}$
bounded random variable with Malliavin derivative given by $D_{s}\xi$
satisfying

$\sup_{0\leq s\leq T}\|D_{s}\xi\|_{L^{2p}}<\infty,\text{ for all }p>1.$

  * (H2)

For all $s\in[0,T]$, $D_{s}\xi$ is Malliavin differentiable, with second order
derivative given by $(D_{s^{\prime}}D_{s}\xi)_{s^{\prime},s\in[0,T]}$ and
satisfying $\sup_{0\leq s^{\prime},s\leq
T}\|D_{s^{\prime}}D_{s}\xi\|_{L^{p}}<\infty$ for all $p>1.$

  * (H3)

$B:\Omega\times[0,T]\rightarrow\mathbb{R}$ and
$C:\Omega\times[0,T]\rightarrow\mathbb{R}^{d}$ are
$\mathfrak{F}_{t}\text{-adapted}$ processes bounded by a constant $\Lambda>0$.
$A:\Omega\times[0,T]\times[0,T]\rightarrow\mathbb{R}$ satisfies $\sup_{0\leq
s\leq T}\|A_{s,\cdot}\|_{\mathcal{S}^{2p}}<\infty$ for all $p>1$. For each
fixed $s\in[0,T]$, the process $(A_{s,t})_{t\in[0,T]}$ is progressively
measurable.

  * (H4)

$A_{s,t},B_{t},C_{t}$ are Malliavin differentiable for all $s\in[0,T]$ and
$t\in[0,T]$. Measurable versions of their Malliavin derivatives are
respectively given by $D_{s^{\prime}}A_{s,t},D_{s}B_{t}$ and $D_{s}C_{t}$ for
$s^{\prime},s\in[0,T]$ such that for all $p>1$

$\sup_{0\leq s^{\prime},s\leq
T}\mathbb{E}\Big{[}\Big{(}\int_{0}^{T}(|D_{s^{\prime}}A_{s,t}|^{2}+|D_{s}B_{t}|^{2}+|D_{s}C_{t}|^{2})\mathrm{d}t\Big{)}^{p}\Big{]}<\infty$

######  Theorem C.1.

Under Assumptions (H1)–(H4), the BSDE (C) has a unique Malliavin
differentiable solution $(U,V)\in\mathcal{S}^{2p}\times\mathcal{H}^{2p}$ for
$p>1$. A version of $(D_{s^{\prime}}U_{s,t},D_{s^{\prime}}V_{s,t})_{0\leq
s^{\prime}\leq s\leq t\leq T}$ satisfies

$\displaystyle D_{s^{\prime}}U_{s,t}$
$\displaystyle=D_{s^{\prime}}D_{s}\xi-\int_{t}^{T}D_{s^{\prime}}V_{s,r}\mathrm{d}W_{r}$
$\displaystyle+\int_{t}^{T}\Big{(}(D_{s^{\prime}}G)(r,s,U_{s,r},V_{s,r})+\langle\nabla
G(r,s,U_{s,r},V_{s,r}),(D_{s^{\prime}}U_{s,r},D_{s^{\prime}}V_{s,r})\rangle\Big{)}\mathrm{d}r.$
(C.2)

Moreover, a version of $(V_{s,t})_{0\leq s\leq t\leq T}$ is given by
$(D_{t}U_{s,t})_{0\leq s\leq t\leq T}$.

Acknowledgments: The authors thank the Editor in Chief and two anonymous
referees for their remarks which helped to considerably improve the
presentation of the paper. A part of this work was carried out while Rhoss
Likibi Pellat visited the Department of Mathematics, Brandenburg Technology
University of Cottbus, whose hospitality is greatly appreciated and he
personally thank Prof Ralf Wunderlich.

## References

  * [1] O. Aboura and S. Bourguin.  Density estimates for solutions to one dimensional backward sde’s.  Potential Analysis, 38:573–587, 2013.
  * [2] F. Antonelli and A. Kohatsu-Higa.  Densities of one-dimensional backward SDEs.  Potential Analysis, 22:263–287, 2005.
  * [3] K. Bahlali.  Solving Unbounded Quadratic BSDEs by a Domination method.  arXiv:1903.11325v1, 03 2019.
  * [4] K. Bahlali, M. Eddahbi, and Y. Ouknine.  Quadratic BSDE with $L^{2}$-terminal data: Krylov’s estimate and Itô-Krylov’s formula and existence results.  Annals of Probability, 45:2377–2397, 2017.
  * [5] V. Bally and L. Caramellino.  Riesz transform and integration by parts formulas for random variables.  Stochastic Processes and their Applications, 121(6):1332–1355, 2011\.
  * [6] D. Banos and P. Krühner.  Hölder continuous densities of solutions of SDEs with measurable and path dependent drift coefficients.  Stoch. Proc. Appli., 88:4:1785–1799, 2017.
  * [7] D. Banos, T. Meyer-Brandis, F. Proske, and S. Duedahl.  Computing deltas without derivatives.  Finance Stoch., 21(2):509–549, 2017.
  * [8] D. Banos and T. Nilssen.  Malliavin and flow regularity of SDEs. Application to the study of densities and stochastic transport equation.  Stochastics, 88:4:540–566, 2016.
  * [9] Y. Chen and P. Luo.  Existence of global solutions for multi-dimensional coupled fbsdes with diagonally quadratic generators.  Journal of Mathematical Analysis and Applications, 522(1):126948, 2023.
  * [10] Z. Chen, R. Kulperger, and G. Wei.  A commonotonic theorem for BSDEs.  Stochastic Process. Appl., 115:41–54, 2005.
  * [11] R. Chertovskih and E. Shamarova.  Gaussian-type Density Bounds for Solutions to Multidimensional Backward SDEs and Application to gene Expression.  Potential Analysis, 2022.
  * [12] F. Delarue.  Estimates of the solutions of a system of quasi-linear PDEs. A probabilistic scheme.  Séminaire de Probabilités XXXVII, 1832:290–332, 2003.
  * [13] F. Delarue and G. Guatteri.  Weak existence and uniqueness for forward-backward sdes.  Stochastic Process. Appl., 116:1712–1742, 2006.
  * [14] J. Detemple, R. Garcia, and M. Rindisbabacher.  Representation formulas for Malliavin derivatives of diffusion processes.  Finance Stoch., 9:349–367, 2005.
  * [15] N.T. Dung.  Gaussian lower bounds for the density via Maliavin calculus.  Comptes Rendus Mathématique, 357:1:79–87, 2020.
  * [16] N.T. Dung, N. Privault, and G.L. Torrisi.  Gaussian estimates for the solutions of some one-dimensional stochastic equations.  Potential Anal., 43:289–311, 2015.
  * [17] N. Eisenbaum.  Integration with respect to local time.  Potential Analysis, 13:303–328, 2000.
  * [18] N. Eisenbaum.  Local time-space stochastic calculus for Lévy processes.  Stoch. Proc. Appl., 116:757–778, 2006.
  * [19] X. Fan and JL. Wu.  Density estimates for the solutions of backward stochastic differential equations driven by gaussian processes.  Potential Analysis, 54:483–501, 2021.
  * [20] F. Flandoli, M. Gubinelli, and E. Priola.  Well posedness of the transport equation by stochastic pertubation.  Invent math, 180:1–53, 2009.
  * [21] F. Flandoli, M. Gubinelli, and E. Priola.  Flow of diffeomorphisms for SDEs with unbounded Hölder continuous drift.  Bull. Sci. math., 134:405–422, 2010.
  * [22] C. Frei and G. Dos Reis.  Quadratic fbsde with generalized burger’s type nonlinearities, pertubations and large deviations.  Stochastics and Dynamics, 13:37 pages, 2013.
  * [23] P. Imkeller, R. Likibi Pellat, and O. Menoukeu Pamen.  Differentiability of quadratic forward-backward SDEs with rough drift.  arXiv:2210.05622v1, 10.
  * [24] P. Imkeller and G. Dos Reis.  Path regularity and explicit convergence rate for bsde with truncated quadratic growth.  Stoch. Process. Appl., 120:348–379, 2010.
  * [25] J. Jackson.  On quasilinear parabolic systems and fbsdes of quadratic growth.  arXiv:2206.08334v2, 04.
  * [26] J.Wei, G. Lv, and W. Wang.  Stochastic transport equation with bounded and dini continuous drift.  J.Diff. Eq., 323, 2022.
  * [27] I. Karatzas and S. E. Shreve.  Brownian Motion and Stochastic Calculus.  Springer-Verlag, New York, 1988.
  * [28] N. Kazamaki.  Continuous exponential martingales and BMO, volume 1579.  Lectures Notes in Mathematics. Springer-Verlag, Berlin.
  * [29] M. Kobylanski.  Backward stochastic differential equations and partial differential equations with quadratic growth.  Annals of Probability, 2:558–602, 2000.
  * [30] N.V. Krylov.  Controlled Diffusion Processes.  Springer-Verlag, New York, 1979.
  * [31] S. Kusuoka.  Existence of densities of solutions of stochastic differential equations by Malliavin calculus.  J. Funct. Anal., 258:758–784, 2010.
  * [32] O.A. Ladyzhenskaya, V.A. Solonnikov, and N.N. Ural’tseva.  Linear and Quasi-linear Equations of Parabolic Type.  American Mathematical Society, Providence, RI, 1968.
  * [33] G. Leoni and M. Morini.  Necessary and sufficient conditions for the chain rule in ${W}_{loc}^{1,1}(\mathbb{R}^{N};\mathbb{R}^{d})$ and ${BV}_{loc}(\mathbb{R}^{N};\mathbb{R}^{d})$.  J. Eur. Math. Soc., 9:219–252, 2005.
  * [34] J. Ma, P. Protter, and J. Yong.  Solving forward-backward stochastic differential equations explicitly: a four step scheme.  Probab. Theory Related Fields, 98:339–359, 1994.
  * [35] J. Ma and J. Yong.  Forward-Backward Stochastic Differential Equations and their Applications.  Springer-Verlag, Berlin, 1999.
  * [36] T. Mastrolia.  Density analysis of non-Markovian BSDEs and applications to biology and finance.  Stoch. Proc. Appl., 128:897–938, 2017.
  * [37] T. Mastrolia, D. Possamai, and A. Reveillac.  Density analysis of BSDEs.  The Annals of Probability, 44:2817–2857, 2016.
  * [38] O. Menoukeu-Pamen, T. Meyer-Brandis, T. Nilssen, F. Proske, and T. Zhang.  A variational approach to the construction and malliavin differentiability of strong solutions of SDE’s.  Math. Ann., 357(2):761–799, 2013.
  * [39] S. Menozzi, A. Pesce, and X. Zhang.  Density and gradient estimates for non degenerate brownian sdes with unbounded measurable drift.  Journal of Differential Equations, 272:330–369, 2021.
  * [40] S. E. A. Mohammed, T. Nilssen, and F. Proske.  Sobolev differentiable stochastic flows for sde’s with singular coefficients: Applications to the stochastic transport equation.  Annals of Probability, 43(3):1535–1576, 2015.
  * [41] D. Nualart.  The Malliavin Calculus and Related Topics.  Springer Berlin, 2nd edition edition, 2006.
  * [42] C. Olivera and E. Shamarova.  Gaussian density estimates for solutions of fully coupled forward-backward of SDEs.  Mathematische Nachtrichten, 293:1554–1564, 2020.
  * [43] C. Olivera and C. Tudor.  Density for solutions to stochastic differential equations with unbounded drift.  Braz. J. Proba. Stat., 33, 2019.
  * [44] O. Menoukeu Pamen and L. Tangpi.  Strong solutions of some one-dimensional SDEs with random. and unbounded drifts.  SIAM J. Appl. Math., 51:4105–4141, 2019.
  * [45] R. Likibi Pellat, O. Menoukeu Pamen, and Y. Ouknine.  A class of quadratic forward-backward stochastic differential equations.  Journal of Mathematical Analysis and Applications, 2022.
  * [46] G. Dos Rei and R.J.N. Dos Rei.  A note on commonotonicity and positivity oof the control components of decoupled quadratic FBSDE.  Stochastic and Dynamics, 13:11(pages), 2013.
  * [47] M. Romito.  A simple method for the existence of a density for stochastic evolutions with rough coefficients.  Electron. J. Probab., 113:1–43, 2018.
  * [48] N. Shintaro.  On pathwise uniqueness and comparison of solutions of one-dimensional stochastic differential equations.  Osaka Journal of Mathematics, 20:197–204, 1983.
  * [49] X. Wang and S. Fan.  A class of stochastic gronwall’s inequality and its application.  Journal of Inequalities and Applications, 2018:336, 2018.
  * [50] A. K. Zvonkin.  A transformation of the state space of a diffusion process that removes the drift.  Math. URSS Sbornik, 22:129–149, 1974.


# First-principles study of the magnetic anisotropy of ultrathin B-, C-, and
N-doped FeCo films

Joanna Marciniak  Mirosław Werwiński  Justyna Rychły-Gruszecka
<EMAIL_ADDRESS>Institute of Molecular Physics,
Polish Academy of Sciences, M. Smoluchowskiego 17, 60-179 Poznań, Poland

###### Abstract

Iron-based layered systems are of great interest because of their ability to
tune effective material parameters such as magnetic anisotropy energy (MAE).
The influence of the crystallographic structure of Fe, its thickness, and the
presence of other layers above and below the Fe layer on magnetic parameters,
such as the MAE of the studied system, is an intriguing and important topic
from an application point of view. Here, we present a density functional
theory (DFT) study of the magnetic anisotropy of nine-monolayer Fe, FeCo, and
FeCo films with B, C, and N dopants placed in octahedral interstitial
positions. The theoretical study is based on calculations using the full-
potential local-orbital code FPLO and the generalized gradient approximation.
The chemical disorder in the FeCo layers was modeled using the virtual crystal
approximation. The structures of the layers were subjected to optimization of
the geometry of the interlayer spacings and the neighborhood of the dopant
sites. We determined the local magnetic moments and the excess charge at each
layer position. We also identified the influence of dopant atoms on the
magnetic properties of FeCo layers, such as magnetization and magnetic
anisotropy.

##  1 Introduction

Permanent magnets are essential in many consumer and industrial products that
convert energy. Rare earth magnets, particularly those based on neodymium and
samarium, typically produce the highest energy product. Of the various
permanent magnets reported to date, NdFeB has the highest magnetic efficiency.
Recently, there has been growing concern about the volatility of the rare
earth market, which manifested itself in the so-called rare earth crisis of
2011 [1]. Over the past decade, the rapidly rising and volatile prices of rare
earth metals have prompted intensive research efforts worldwide to develop
alternative materials for permanent magnets, especially those with very low or
zero rare earth content. FeCo has the highest saturation magnetization of all
transition metals and their alloys. Still, due to its cubic crystal structure
(bcc), it has an extremely low magnetic anisotropy, making it a typical soft
magnetic material, which has long been unsuitable for permanent magnets.
Burkert et al. [2] performed first-principles calculations using a virtual
atom model to theoretically predict that FeCo with a tetragonal crystal
structure (bct) would exhibit a giant magnetic anisotropy constant of up to
$10^{7}~{}\mathrm{J\,m}^{-3}$. In a FeCo alloy with lattice parameter ratio
$c/a\approx 1.2$ indication of high tetragonal deformation, the magnetic
anisotropy energy (MAE) was predicted to reach $0.7-0.8$ meV atom-1. However,
the experimental growth of Fe-Co alloys with such high tetragonal distortion
is challenging because their natural structure is bcc (up to about 70% Co
concentration). Subsequent theoretical studies have corrected the theoretical
MAE to values very close to the experimental ones, using a more realistic
treatment of the chemical disorder present in Fe-Co alloys [3, 4]. To solve
this problem, tiny interstitial atoms (B, C, N, etc.) placed on the $c$-axis
of FeCo can stretch the unit cell and stabilize the bct structure. As for
FeCo, there are studies on adding interstitial elements B [5, 6, 7], C [8, 9,
10], N [11, 12, 13].

In the race to miniaturize electronic devices, modern science has reached the
ultimate level of monoatomic layers. The limit for magnetic tunnel junctions
with perpendicular magnetic anisotropy is about twenty atomic monolayers [14,
15]. Layered systems are particularly interesting for their ability to tune
effective material parameters such as MAE. Among such systems, iron-based
layered systems are of considerable interest [16, 17, 18, 19]. An intriguing
and important topic from the point of view of applications is the influence of
the crystallographic structure of Fe, its thickness, and the presence of other
layers above and below the Fe layer on magnetic parameters such as the MAE of
the studied system.

Figure 1:  A unit cell of a 9-atomic monolayer thick Fe0.7Co0.3 thin film,
with the dopant (B or C or N) located in an octahedral interstitial position.
A single computational cell is shown on the left, and a multiplied $2\times 2$
cell on the right. The unit cell is obtained after geometry optimization. The
red spheres represent the Fe0.7Co0.3 virtual atoms, while the green sphere
represents the dopant atom. A surrounding vacuum (52 Å) was applied in the
computational cells to obtain non-interacting thin films. To improve the
visibility of the film’s atoms, the height of the vacuum used has been
significantly reduced in the figures.

Here, we present a theoretical study of the magnetic anisotropy of FeCo thin
films with B, C, and N dopants located in octahedral interstitial positions.
The theoretical study is based on calculations using the full-potential local-
orbital electronic structure code (FPLO) and the generalized gradient
approximation. The chemical disorder in FeCo layers is modeled using the
virtual crystal approximation (VCA). The layer structures are subjected to
geometry optimization of the interlayer distances and the vicinity region of
the dopant sites. We determine the local magnetic moments and excess charge at
each atomic position. We identified the dopant atoms’ effect on the FeCo
films’ magnetic properties, such as magnetization and magnetic anisotropy.

##  2 Calculations’ details

Figure 2: Magnetic anisotropy energy (MAE) calculated for 9 atomic monolayer
thick films. From left to right: pure Fe (shown as a purple bar); pure
Fe0.7Co0.3 alloy (shown as a green bar); Fe0.7Co0.3 with a B dopant atom
(shown as a red bar); Fe0.7Co0.3 with a C dopant atom (shown as a blue bar);
Fe0.7Co0.3 with an N dopant atom (shown as a yellow bar). The dopant atoms are
in the octahedral interstitial position. The DFT calculations were performed
with the full-potential local-orbital (FPLO-18) code [20]. The exchange-
correlation functional in the Perdew-Burke-Ernzerhof (PBE) parametrization was
used [21].

In this study, we examined the structural and magnetic properties of ultrathin
(9 atomic monolayers thick) films of Fe0.7Co0.3 alloys with a dopant (of B, C,
and N atoms) located in octahedral interstitial position, see Fig. 1, using
density functional theory (DFT) calculations. To compare the effect of the
dopant, placed in the octahedral gap, on the obtained magnetic properties of
the structures, we also performed calculations of thin films of pure iron-
cobalt alloy with a given concentration (70% Fe and 30% Co) and a thin film of
iron without doping. Altogether, calculations of five different structures
were performed.

The structures were prepared in cells with symmetry $P4/mmm$, where the
lattice constant equals $a=4.0~{}\AA{}$. In the computational cells, to remove
the bulk periodicity in the indicated direction and to model the thin film, a
vacuum of 52 Å was added, giving a height of the computational cell equal to
64 Å. In the presented calculations, the VCA approach was used, as implemented
in FPLO 18.00-52 code, to simulate a disordered 30% Co alloy. Wyckoff’s
positions were optimized in each considered system using a scalar-relativistic
approach with spin polarization. Atomic positions in [001] were optimized by
forces optimization on atoms, with convergence criterion set as 10-3 eV Å-1.

A mesh of $k$-points was set as 65 $\times$ 65 $\times$ 5, and the convergence
criteria were used 10-8 Ha and 10-7 for energy and charge density,
respectively. The integration over the Brillouin zone was performed using the
tetrahedron method. The DFT calculations were performed using the full-
potential local-orbital approach implemented in the FPLO 18.00-52 code [20].
Using the Perdew-Burke-Ernzerhoff (PBE) exchange-correlation potential [21],
scalar-relativistic calculations were performed, followed by 1 iteration of
full-relativistic calculations. MAE was determined by the formula:

$\mathrm{MAE}=E(\theta=0^{\circ})-E(\theta=90^{\circ}),$ (1)

where $\theta$ is the angle between the magnetization direction and the $c$
axis. Crystal structures were visualized using the VESTA program [22].

##  3 Results and Discussion

Table 1:  Total thicknesses of the 9-atomic monolayer thick Fe0.7Co0.3 and
Fe0.7Co0.3-X thin films after interlayer spacing optimization, together with
distances between the central atomic monolayer and the first and second-
nearest monolayer. The double values for the atoms of the first nearest
monolayer are due to the presence of dopant atoms in the central layer (see
Fig. 1). The $a$ lattice parameter of the computational unit cell was set at
4.00 Å ($\sqrt{2}\times 2.83$ Å).
| (Fe0.7Co0.3)18 | (Fe0.7Co0.3)18B | (Fe0.7Co0.3)18C | (Fe0.7Co0.3)18N
---|---|---|---|---
total thickness (Å) | 11.27 | 11.81 | 11.74 | 11.75
1st layer (Å) | 1.42 | 1.58/1.86 | 1.59/1.80 | 1.60/1.80
2nd layer (Å) | 2.83 | 3.09 | 3.06 | 3.07
Figure 3: Excess charge (charge transfer) at each atomic position in the calculated thin films presented as black empty squares and big black full triangles for dopant atoms. Local spin magnetic moments at each atomic position in the calculated thin films presented as red empty diamonds and big red full circles for dopant atoms. Considered 9-atomic monolayer thick films of (a) pure Fe; (b) pure Fe0.7Co0.3 alloy; (c) Fe0.7Co0.3 with a B dopant atom; (d) Fe0.7Co0.3 with a C dopant atom; (e) Fe0.7Co0.3 with an N dopant atom. The dopant atoms in (c-e) are in the octahedral interstitial position. Calculations were performed with the PBE exchange-correlation potential [21] in scalar-relativistic formalism using the FPLO18 code [20].  Figure 4:  The calculated (a) spin and (b) orbital magnetic moment for all considered 9-atom layers (from left to right: pure Fe; pure Fe0.7Co0.3 alloy; Fe0.7Co0.3 with B-, C-, or N-atom, located in the octahedral interstitial position, respectively. The red, yellow, and orange bars represent the (a) spin (b) orbital magnetic moment calculated for the outer, central, average atom (Fe in case of pure Fe thin film or Fe0.7Co0.3 for the rest of the structures), respectively; the green bar represents a (a) spin (b) orbital magnetic moment for the dopant atom. The DFT calculations were performed with the full-potential local-orbital (FPLO-18) code [20]. The exchange-correlation functional in the Perdew-Burke-Ernzerhof (PBE) parametrization was used [21]).  Figure 5:  The calculated total density of states (DOS) for the considered 9-atom thick layers of: (a) pure Fe; (b) pure Fe0.7Co0.3 alloy; (c) Fe0.7Co0.3 with a B dopant atom; (d) Fe0.7Co0.3 with a C dopant atom; (e) Fe0.7Co0.3 with an N dopant atom. The dopant atoms in (c-e) are in the octahedral interstitial position. The DFT calculations were performed with the full-potential local-orbital (FPLO-18) code [20]. The exchange-correlation functional in the Perdew-Burke-Ernzerhof (PBE) parametrization was used [21].  Table 2: Averaged spin ($m_{s}$) and orbital ($m_{l}$) magnetic moments, together with magnetocrystalline anisotropy energies (MAE) of 9-monolayer thick films of Fe, Fe0.7Co0.3, and Fe0.7Co0.3X models, where X = B, C, N. The obtained results were compared with literature values for similar systems.  System | Form | Method | Reference | $m_{s}$ | $m_{l}$ | MAE | MAE
---|---|---|---|---|---|---|---
|  |  |  | ($\mu_{\mathrm{B}}$ at-1) | ($\mu_{\mathrm{B}}$ at-1) | ($\mu$eV at-1) | (MJ m-3)
Fe | 9-monolayers |  | this work | 2.41 | 0.058 | 87 | 1.38
Fe | bulk |  | [15] | 2.16 | 0.043 | 0 | 0
Fe0.7Co0.3 | 9-monolayers | VCA | this work | 2.38 | 0.080 | 43 | 0.68
Fe0.7Co0.3 | bulk | VCA/CPA | [9] | 2.42/2.43 | – | 54/27 | –
(Fe0.7Co0.3)18B | 9-monolayers | VCA | this work | 2.18 | 0.076 | -13 | -0.18
(Fe0.65Co0.35)24B | bulk | CPA | [5] | 2.16 | 0.064 | 23 | 0.31
Fe67Co18B15 | >100 nm film | expt. | [7] | $\sim 1.7$ | – | $\sim 0$ | $\sim 0$
(Fe0.38Co0.62)0.98B0.02 | 20 nm films | expt. | [5] | 1.9 $\pm$ 0.2 | – | – | 0.4 $\pm$ 0.2
(Fe0.7Co0.3)18C | 9-monolayers | VCA | this work | 2.16 | 0.075 | 6 | 0.09
(Fe0.7Co0.3)8C | bulk | VCA/CPA | [9] | 1.88/1.85 | – | 77/81 | 1.10/–
(Fe0.69Co0.31)16C | bulk | supercells | [10] | 2.35 | – | 65 | 0.95
(Fe0.4Co0.6)0.98C0.02 | 5 nm films | expt. | [8] | 2.0 | – | – | 0.80 $\pm$ 0.15
(Fe0.7Co0.3)18N | 9-monolayers | VCA | this work | 2.19 | 0.076 | 22 | 0.33
(Fe0.75Co0.25)16N2 | bulk | supercells | [13] | 1.96 | – | – | 1.5
(Fe0.5Co0.5)89V9N2 | 20 nm film | expt. | [11] | – | – | – | 1.24

We started the calculations for the considered structures by optimizing the
$z$ positions of all atomic sites with a fixed value of the lattice parameter
$a$. The total thicknesses of the Fe-Co thin films and the distances of the
nearest and second-nearest atomic monolayers to the central layer were
collected in Table 1. As we can see, the film’s total thickness increases due
to doping. The lattice stretching in the $z$ direction is caused by the dopant
atom located in the octahedral hole.

After that, we calculated the magnetic anisotropy energy of FeCo thin films (9
atoms thick) with B, C, and N dopants in octahedral interstitial position in
the center of the layer. For comparison, we have calculated a pure thin film
of an iron-cobalt alloy of a given concentration and a thin film of iron. The
MAE for the above structures is shown in Fig. 2. For the structures of a given
thickness, we can see that pure Fe layers have the strongest tendency to be
magnetized perpendicular to the surface; in the case of the Fe0.7Co0.3 alloy,
the MAE decreases almost by half, while for the thin films doped with N, C,
and B elements the MAE becomes smaller and smaller. In the case of the
Fe0.7Co0.3 with a B dopant atom, the MAE changes sign, which shows the change
in the tendency of this structure and its preference to be magnetized in-plane
(and not out-of-plane as was the case for the other systems).

While in-plane anisotropy resulting from B doping appears uncommon in
comparison to other cases, it is frequently obtained in FeCoB thin films
through experimental means [23]. Based on earlier experimental work, the
orientation of magnetization in FeCoB films can be altered through various
technological methods, including post-deposition annealing in a magnetic field
or application of mechanical stress [23]. The direction of magnetization is
influenced by various factors, such as the concentration of individual
elements, the type and thickness of the substrate and overlay, and the
thickness of the magnetic film. The calculated result for the ultrathin
Fe0.7Co0.3-B film implies that it has an intrinsic property of in-plane
magnetic anisotropy, unaffected by the aforementioned factors.

We also have conducted a detailed study of the charge and the spin magnetic
moment for all the systems under study and present the results in Fig. 3. We
first plotted (for comparison) the charge transfer at each atom of pure Fe and
pure Fe0.7Co0.3 alloy through the thickness of their layers using black empty
squares connected by a black line (see Fig. 3(a) and (b), respectively). Only
negligibly small charge oscillations can be observed for these single-
component layers, appearing due to the presence of surfaces. By introducing a
dopant into the octahedral gap located in the center of a thin FeCo alloy
layer, we can observe the transfer of a negative charge by the dopant atom at
the expense of the alloy atoms in its vicinity (see Fig. 3(c-e) for dopants
with B, C, and N, respectively, where the additional charge on the dopant has
been visualized with a full black triangle).

We then plotted, using red empty diamonds connected by a red line, the local
spin magnetic moments at each atomic position in the calculated thin films
(see Fig. 3). It can be observed that the introduction of surfaces into the
structure of pure magnetic layers without dopants increases the spin magnetic
moment near their surface. This is particularly evident for a thin layer of
pure Fe (see Fig. 3(a)), a slightly less pronounced increase in the magnetic
moment is obtained for a thin layer of Fe0.7Co0.3 alloy (see Fig. 3(b)). For
comparison, the thin dashed line has been used to indicate the values of the
spin magnetic moment for solid Fe
($m_{\rm{s}}~{}=~{}2.2~{}\mu_{\rm{B}}\rm~{}{at}^{-1}$, marked in Fig. 3(a))
and for solid Fe0.7Co0.3
($m_{\rm{s}}~{}=~{}2.42~{}\mu_{\rm{B}}\rm~{}{at}^{-1}$, marked in Fig.
3(b-e)). By introducing the dopant into the octahedral gap in the center of
the Fe0.7Co0.3 alloy thin film, we can observe a varying spin magnetic moment
on the alloy atoms near the dopant. The dopant atom acquires a small negative
spin magnetic moment (see Fig. 3(c-e) for dopants with B, C, and N,
respectively, where the spin magnetic moment on the dopant has been visualized
with a large red solid circle). The values of spin and orbital magnetic
moments on the dopants obtained in the calculations are summarized in Table 3.
Averaged values of magnetic moments at Fe/Co position and magnetic moments
from literature are also summarized in Table 2. As can be seen, our results
are comparable to those obtained in experiments on thin films and calculations
for solids.

To make the spin magnetic moment results easier to read, we have presented the
obtained results in the form of a bar chart, where the averaged spin magnetic
moment is shown with an orange bar, the spin magnetic moment for the central
(outer) atom is shown with a yellow (red) bar, and the spin magnetic moment
for the dopant atom is shown with a green bar (see Fig. 4(a)). Similarly, we
presented the results for the orbital magnetic moment (see Fig. 4(b)).

Table 3:  Spin and orbital magnetic moments (in $\mu_{B}$) on X = B, C, and N in researched systems. The DFT calculations were performed with the full-potential local-orbital (FPLO-18) code [20]. The exchange-correlation functional in the Perdew-Burke-Ernzerhof (PBE) parametrization was used [21].  System | $m_{(\rm{s},X)}$ | $m_{(\rm{l},X)}$
---|---|---
(Fe0.7Co0.3)18B | -0.33 | 0.002
(Fe0.7Co0.3)18C | -0.25 | 0.002
(Fe0.7Co0.3)18N | -0.03 | 0.004

At the end, we present the calculation of the total density of states to show
what happens after doping (Fe0.7Co0.3)18 thin films with B, C or N, see Fig.
5. First, in Fig. 5(a), we show the total DOS for the pure Fe thin film with
the comparison to the pure (Fe0.7Co0.3)18 thin film shown in Fig. 5(b).
Compared to pure Fe, (Fe0.7Co0.3)18 has 0.3 more electrons per transition
metal atom than pure Fe, causing the Fermi level to shift slightly to the
right for the (Fe0.7Co0.3)18 alloy. Doping the (Fe0.7Co0.3)18 alloy layer with
B, C, or N atoms (see Fig. 5(c-e)) results in a small additional contribution
of $p$-group electrons, which causes a further shift of the Fermi level to the
right. This shift is slight because there is one dopant atom for every 18
atoms of $4s$ and $3d$ band atoms in the (Fe0.7Co0.3)18 alloy thin film. It
can be seen that the dopants form a bond with the alloy, as we do not observe
sharp peaks in the total DOS curve.

Finally, we would like to comment on the selected limitations of the discussed
models resulting from the choice of the size of the system and the usage of
virtual crystal approximation. The thickness of the 9-monolayers was chosen
slightly above the experimentally observed structural transition of ultrathin
Fe films from a structure with strong tetragonal deformation ($c/a$ > 1.4)
below the critical thickness to a bcc-like structure above it [24]. At the
same time, we did not choose a higher thickness to not increase the number of
non-equilibrium atoms in the model and thus the time required for
calculations. Since we did not consider a range of film thicknesses (e.g.,
from 1 to 20 atomic monolayers), the value of our results for a single unique
thickness is limited. For example, in line with the authors’ other research on
ultra-thin films [25, 26] the variation in MAE with thickness is not monotonic
and can change significantly with a change in thickness by as little as one or
two monolayers. Hence, the most significant conclusion regarding the MAE of
the examined layers is the tendency of MAE to increase with the increase in
the atomic number of the dopant. On the other hand, the calculated magnetic
moments on the film’s surface, around the dopant, and on it depend not so much
on the film thickness.

The prepared models also have a specific location of the dopant to the other
atoms - it is part of the central layer of the system, creating a very regular
distribution of the dopant in it, see Fig. 1. This assumption was dictated by
the need to limit the unit cell size, counted in the number of non-equivalent
atoms, to a value enabling MAE calculations. However, placing dopant atoms in
a single monolayer primarily results in overestimating the increase of film
thickness due to doping. However, as experiments on FeCo thin films with B and
C dopants have shown, the magnitude of the tetragonal deformation in our
models is relatively consistent with the measurements [8, 5].

Another approximation used in the presented study is the virtual crystal
approximation (VCA), which we applied to model Fe/Co alloying. We decided to
use it even though it is known for overestimating the value of the MAE for
bulk materials. In the case of modeling layered systems containing
interstitial defects, which by breaking the symmetry leads to a significant
increase in the number of non-equivalent atoms in the system and thus to a
significant increase in computational time, the use of the VCA allows both to
optimize the system geometry and to reduce the computational time by at least
an order of magnitude, relative to alternative methods of (1) averaging over
multiple supercells or (2) the coherent potential approximation (CPA), the
latter usually happens in combination with the atomic sphere approximation
(ASA), leading to unreliable MAE results.

Comparing the obtained MAE values with the results from the literature for
calculations for bulk systems (see Table 2), we notice that the obtained
values differ significantly. However, these differences may primarily result
from the fact that calculations of the layered system are performed, where
size effects come into play. Even calculations provided using the VCA method
in the bulk system resulted in much higher values. Comparison with the results
of experiments on thin-film systems allows us to confirm the similarity of
orders of magnitude.

##  4 Summary and Conclusions

We have conducted a theoretical study of the magnetic anisotropy of ultrathin
(9 atomic monolayers thick) FeCo films with B, C, and N dopants located in
octahedral interstitial position in the center of the layer. The layer
structures were subjected to geometry optimization of the interlayer distances
and the vicinity region of the dopant sites. This allowed us to get different
distances between different layers according to the placement above/below the
dopant atom or next to this place and in the different depths of the film
thickness. We determined the local magnetic moments and excess charge at each
position in the films and identified the dopant atoms’ effect on the FeCo
films’ magnetic properties, such as magnetization, total DOS, and magnetic
anisotropy.

Contrary to the results for bulk systems, our results indicate that doping
FeCo with B, C, and N atoms in the octahedral position in 9 atomic monolayers
thick FeCo films can significantly reduce the MAE, even changing its sign to
negative in the case of B-doped FeCo thin films. The results of this
investigation may have important implications for further research on magnetic
thin films for spintronic applications.

## Acknowledgements

We acknowledge the financial support of the National Science Centre Poland
under the decision DEC-2018/30/E/ST3/00267. We thank Paweł Leśniak and Daniel
Depcik for compiling the scientific software and administration of the
computing cluster at the Institute of Molecular Physics, Polish Academy of
Sciences.

## References

  * [1] J. Coey, Permanent magnets: Plugging the gap, Scr. Mater. 67 (6) (2012) 524–529.  doi:https://doi.org/10.1016/j.scriptamat.2012.04.036.
  * [2] T. Burkert, L. Nordström, O. Eriksson, O. Heinonen, Giant magnetic anisotropy in tetragonal FeCo alloys, Phys. Rev. Lett. 93 (2) (2004) 027203.  doi:10.1103/PhysRevLett.93.027203.
  * [3] C. Neise, S. Schönecker, M. Richter, K. Koepernik, H. Eschrig, The effect of chemical disorder on the magnetic anisotropy of strained Fe–Co films, Phys. Status Solidi (b) 248 (10) (2011) 2398–2403.  doi:10.1002/pssb.201147100.
  * [4] I. Turek, J. Kudrnovskỳ, K. Carva, Magnetic anisotropy energy of disordered tetragonal Fe-Co systems from ab initio alloy theory, Phys. Rev. B 86 (17) (2012) 174430.  doi:10.1103/PhysRevB.86.174430.
  * [5] L. Reichel, L. Schultz, D. Pohl, S. Oswald, S. Fähler, M. Werwiński, A. Edström, E. K. Delczeg-Czirjak, J. Rusz, From soft to hard magnetic Fe–Co–B by spontaneous strain: a combined first principles and thin film study, J. Phys. Condens. Matter 27 (47) (2015) 476002.  doi:10.1088/0953-8984/27/47/476002.
  * [6] W. Coene, F. Hakkens, R. Coehoorn, D. De Mooij, C. De Waard, J. Fidler, R. Grössinger, Magnetocrystalline anisotropy of Fe3B, Fe2B and Fe1.4Co0.6B as studied by Lorentz electron microscopy, singular point detection and magnetization measurements, J. Magn. Magn. Mater. 96 (1-3) (1991) 189–196.  doi:10.1016/0304-8853(91)90627-M.
  * [7] L. Chen, T. Klemmer, K. Ellis, R. Van Dover, S. Jin, Soft-magnetic properties of Fe–Co–B thin films for ultra-high-frequency applications, J. Appl. Phys. 87 (9) (2000) 5858–5860.  doi:10.1063/1.372546.
  * [8] L. Reichel, L. Schultz, S. Fähler, Lattice relaxation studies in strained epitaxial Fe-Co-C films, J. Appl. Phys. 117 (17) (2015).  doi:10.1063/1.4908031.
  * [9] E. K. Delczeg-Czirjak, A. Edström, M. Werwiński, J. Rusz, N. V. Skorodumova, L. Vitos, O. Eriksson, Stabilization of the tetragonal distortion of FexCo1-x alloys by C impurities: A potential new permanent magnet, Phys. Rev. B 89 (2014) 144403.  doi:10.1103/PhysRevB.89.144403.
  * [10] W. Marciniak, M. Werwiński, Structural and magnetic properties of Fe-Co-C alloys with tetragonal deformation: a first-principle study (2023).  arXiv:2307.05709, doi:10.48550/arXiv.2307.05709.
  * [11] T. Hasegawa, T. Niibori, Y. Takemasa, M. Oikawa, Stabilisation of tetragonal FeCo structure with high magnetic anisotropy by the addition of V and N elements, Sci. Rep. 9 (1) (2019) 5248.  doi:10.1038/s41598-019-41825-7.
  * [12] T. Hasegawa, Challenges toward development of rear-earth free FeCo based permanent magnet, Electron. Commun. Jpn 104 (2) (2021) e12307.  doi:10.1002/ecj.12307.
  * [13] X. Zhao, C.-Z. Wang, Y. Yao, K.-M. Ho, Large magnetic anisotropy predicted for rare-earth-free Fe16-xCoxN2 alloys, Phys. Rev. B 94 (22) (2016) 224424.  doi:10.1103/PhysRevB.94.224424.
  * [14] K. Watanabe, B. Jinnai, S. Fukami, H. Sato, H. Ohno, Shape anisotropy revisited in single-digit nanometer magnetic tunnel junctions, Nat. Commun. 9 (1) (2018) 663.  doi:10.1038/s41467-018-03003-7.
  * [15] J. Snarski-Adamski, J. Rychły, M. Werwiński, Magnetic properties of 3d, 4d, and 5d transition-metal atomic monolayers in Fe/TM/Fe sandwiches: Systematic first-principles study, J. Magn. Magn. Mater. 546 (2022) 168828.  doi:10.1016/j.jmmm.2021.168828.
  * [16] A. Perumal, Y. Takahashi, T. Seki, K. Hono, Particulate structure of L10 ordered ultrathin fept films for perpendicular recording, Appl. Phys. Lett. 92 (13) (2008).  doi:10.1063/1.2830708.
  * [17] K. Hammar, Y. Labaye, L. Messad, A. Ziane, Theoretical estimation of surface magnetic anisotropy on L10-FePt thin films: Case of perfect and defect surfaces., Surf. Sci. 717 (2022) 121999.  doi:10.1016/j.susc.2021.121999.
  * [18] Y. Miura, S. Ozaki, Y. Kuwahara, M. Tsujikawa, K. Abe, M. Shirai, The origin of perpendicular magneto-crystalline anisotropy in L10–FeNi under tetragonal distortion, J. Phys. Condens. Matter 25 (10) (2013) 106005.  doi:10.1088/0953-8984/25/10/106005.
  * [19] K. Ito, M. Hayashida, H. Masuda, T. Nishio, S. Goto, H. Kura, T. Koganezawa, M. Mizuguchi, Y. Shimada, T. J. Konno, H. Yanagihara, K. Takanashi, Epitaxial L10–FeNi films with high degree of order and large uniaxial magnetic anisotropy fabricated by denitriding FeNiN films, Appl. Phys. Lett. 116 (24) (2020).  doi:10.1063/5.0011875.
  * [20] K. Koepernik, H. Eschrig, Full-potential nonorthogonal local-orbital minimum-basis band-structure scheme, Phys. Rev. B 59 (3) (1999) 1743–1757.  doi:10.1103/PhysRevB.59.1743.
  * [21] J. P. Perdew, K. Burke, M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77 (18) (1996) 3865–3868.  doi:10.1103/PhysRevLett.77.3865.
  * [22] K. Momma, F. Izumi, _VESTA_ : A three-dimensional visualization system for electronic and structural analysis, J. Appl. Crystallogr. 41 (3) (2008) 653–658.  doi:10.1107/S0021889808012016.
  * [23] L. Kipgen, H. Fulara, M. Raju, S. Chaudhary, In-plane magnetic anisotropy and coercive field dependence upon thickness of CoFeB, J. Magn. Magn. Mater. 324 (19) (2012) 3118–3121.  doi:10.1016/j.jmmm.2012.05.012.
  * [24] B. R. Cuenya, M. Doi, S. Löbus, R. Courths, W. Keune, Observation of the fcc-to-bcc Bain transformation in epitaxial Fe ultrathin films on Cu3Au (001), Surf. Sci. 493 (1-3) (2001) 338–360.  doi:10.1016/S0039-6028(01)01239-0.
  * [25] J. Marciniak, M. Werwiński, L10 FePt thin films with tilted and in-plane magnetic anisotropy: First-principles study (2023).  arXiv:2303.15877, doi:10.48550/arXiv.2303.15877.
  * [26] J. Meixner, J. Rychły-Gruszecka, M. Werwiński, Magnetic properties and structural phase transition in ultrathin fcc Fe(111) and bccFe (111) films: first-principles study, arXiv preprint arXiv:2307.11167 (2023).  doi:10.48550/arXiv.2307.11167.