\section{Introduction}
\label{sec:intro}
Inside a cell, actin filaments grow and exert polymerization force on the plasma membrane and create membrane protrusions. This process plays an important role in cell motility \cite{review1, review2, review3, review4, review5,farrel2013,banguy2013}. In many experiments, it has been studied how the mechanical properties of the plasma membrane directly affect the formation of membrane protrusion. For example, by artificially increasing the membrane tension in experiments, it has been seen that the rate of protrusion formation goes down, while by decreasing the membrane tension, the rate is found to increase \cite{NilsGauthier2011, DRaucher2000}. Although, these studies indicate that the membrane tension is generally an obstacle to protrusion formation, in \cite{EllenLBatchelder2011}, it was shown that in certain situations, membrane tension can also enhance protrusions by streamlining actin polymerization in one specific direction.

These studies show it is important to understand how the elastic interactions in the membrane influence the mechanism of protrusion formation and in this rapid communication, we address this question in a simple setting. We model the elastic membrane using the Helfrich Hamiltonian, which is the most commonly used model that includes surface and bending energy of the membrane \cite{Helfrich1973, Julian2016, Lawrence2005, Daillant2005, Quemeneur2014, Bossa2018, Fournier2014, Lenz2009, lin2006, peleg2011, mark2010, isaac2013, shlomovitz2011, veksler2007, veksler2009, orly2014}. We describe the membrane using a one dimensional height field, whose time-evolution is governed by the Helfrich Hamiltonian. Any variation in the height costs energy and a flat membrane corresponds to the lowest energy configuration. We consider the membrane being pushed by few growing filaments, which tend to create protrusions in the membrane that cost energy \cite{baumgaertner2010, baumgaertner2012, Sadhu2018, lipowsky94}. As surface tension $\sigma$ or bending rigidity $\kappa$ of the membrane is increased, one would expect that protrusion formation should become more difficult, since the energy cost for creating a protrusion ought to go up monotonically with $\sigma$ and $\kappa$. Surprisingly, we find it is not so. For a fixed value of $\kappa$, there is a significant range of $\sigma$, for which the energy cost actually decreases with $\sigma$, which makes protrusion formation easier. More specifically, the surface energy cost for creating a local protrusion in the membrane increases monotonically with $\sigma$, as expected, but the bending energy cost shows a peak, which in turn gives rise to a peak in the total energy cost. Not only for cell motility, our finding has potential implications for a wide class of systems, where elastic deformation of a membrane is involved. The fact that we have been able to observe this effect in a simple and general model is encouraging and this opens up the exciting possibility of finding it in many different kinds of systems.

To understand the mechanism behind this intriguing effect, we examine the shape of the membrane near the binding site, where the filament is in contact with the membrane, and show that a qualitative change in the shape is responsible for this. In Fig. \ref{avg_height_schematic}, we depict this mechanism. In the limit when both $\sigma$ and $\kappa$ are large, the membrane has very slow spatial variation of its height. In this limit, one can neglect higher order derivatives of height and assume the height profile of a local protrusion is almost linear in space, as shown in Fig. \ref{avg_height_schematic} by the thin red line. For such a local shape of the membrane, the bending energy (which scales as the square of the second derivative of height) has a non-zero value at the binding site. Everywhere else on the membrane near the binding site, the bending energy is negligible. As $\sigma$ decreases slightly but still remains large, the shape of the membrane remains qualitatively the same, but the height gradient magnitude now increases (see the blue dotted line in Fig. \ref{avg_height_schematic}) and the bending energy becomes significantly larger at the binding site. To minimize this energy, as $\sigma$ is lowered further, the membrane shape finally changes qualitatively, and becomes as shown by the black thick line in Fig. \ref{avg_height_schematic}, where the peak gets rounded and the height profile also does not remain linear anymore. While the bending energy now is non-zero even away from the binding site, its variation across the membrane happens more gradually. The bending energy cost to create a local protrusion, which is proportional to the fourth derivative of height, is lower for such a configuration.
\begin{figure}[h!]
\centering
\includegraphics[scale=0.7]{avg_height_scematic.eps}
\caption{Schematic picture of typical membrane shape $h(x)$ around the binding site at $x=b$. Keeping $\kappa$ constant at a moderate or large value, as $\sigma$ is varied, the membrane shape changes qualitatively. The curve at the bottom (top) is for the largest (smallest) $\sigma$ value.}
\label{avg_height_schematic} \end{figure}

\section{Description of the model}
\label{sec:model}
Our system consists of a set of $N$  parallel filaments growing against a membrane as shown in Fig. \ref{model}. The membrane is described by a height profile $\{h_i\}$, defined on a one dimensional lattice of length $L$ with lattice constant unity. We assume periodic boundary condition $h_i = h_{i+L}$ on the lattice. In the absence of any external force, the height profile of the membrane follows an equilibrium dynamics with the Helfrich Hamiltonian, which is standardly used to describe the height fluctuations of a plasma membrane \cite{lin2006,mark2010,isaac2013,shlomovitz2011}. This Hamiltonian consists of surface interaction and bending interaction of the membrane and in our lattice model, it has the form \cite{nelson, lipowsky94, AVolmer1998, baumgaertner2012}
\begin{equation}
{\cal H}=\sigma \sum_{i=1}^L (h_i')^2 + \kappa \sum_{i=1}^L (h_i'')^2=\sigma \sum_{i=1}^L (h_i-h_{i-1})^2 + \kappa \sum_{i=1}^L (h_{i-1}-2h_i+h_{i+1})^2,
\label{eq:hamiltonian}
\end{equation}
where $\sigma$ is proportional to the surface tension and $\kappa$ to the bending rigidity. Both of these parameters have dimensions of energy/length$^2$, as follows from the above equation. In Eq. (\ref{eq:hamiltonian}), we have neglected the non-linear terms, which can be justified if the magnitude of height gradient everywhere on the membrane is much less than unity \cite{veksler2007,nirgov2006,Kabaso2011}. The size of a typical eukaryotic cell is about a few tens of micrometers \cite{book2}, which can be compared with the size of our membrane patch. The horizontal distance between two consecutive lattice sites in our model can be taken to be $1/L$ times the size of the membrane. The vertical height difference between two consecutive lattice sites is in the scale of nanometers, set by the size of an actin monomer. Thus, height gradient remains small and the linear approximation for the Hamiltonian in Eq. (\ref{eq:hamiltonian}) remains valid.

It follows from Eq. (\ref{eq:hamiltonian}) that the minimum energy configuration is reached when the membrane is completely flat, {\sl i.e.} $h_i$ the same for all $i$. For a finite temperature, the membrane undergoes thermal fluctuations, as a result of which $h_i$ can increase or decrease by an amount $\delta$. The corresponding energy costs $\mathcal{E}_i ^\pm$ can be obtained from  Eq. (\ref{eq:hamiltonian}) as
$$\mathcal{E}_i^\pm=\sigma \{2 \delta ^2 \mp 2 \delta (h_{i-1}-2h_i+h_{i+1})\} + \kappa \{6 \delta ^2 \pm 2 \delta (h_{i-2}-4h_{i-1} +6h_i-4h_{i+1} + h_{i+2})\}$$
or, in terms of the discrete derivatives,
\begin{equation}
\mathcal{E}_i^\pm=\sigma \{2 \delta ^2 \mp 2 \delta h_i''\} + \kappa \{6 \delta^2 \pm 2 \delta h_i''''\}= \Sigma_i^ \pm  + \mathcal{K}_i^\pm.
\label{eq:energy_cost} \end{equation}
Here $\mathcal{E}_i ^ +$ denotes the total energy cost for changing $h_i$ to $h_i + \delta$ and $\Sigma_i^ +$ denotes the surface energy cost and $\mathcal{K}_i^+$ denotes the bending energy cost for the same process. Similarly,  $\mathcal{E}_i ^ -$, $\Sigma_i^ -$, and $\mathcal{K}_i^-$  are the respective energy costs for changing $h_i$ to $h_i -\delta$. We use the Metropolis algorithm to perform the simulations.

One additional constraint that the system must obey is that the height of the membrane at the binding sites should be such that the membrane always stays above the filament tips (as shown in Fig. \ref{model}). Any height fluctuation that brings a binding site at a lower height than the filament tip, is forbidden. As long as this constraint is satisfied, the update rules can be chosen following local detailed balance $R_+ / R_- = e^{-\beta \epsilon}$, where $R_+$ is the rate of a process that has a positive energy cost $\epsilon$ and $R_-$ is the rate of the reverse process.

\begin{figure}[h!]
\centering
\includegraphics[scale=0.8]{model_bending.eps}
\caption{Schematic diagram of the model. The thick solid line represents the height profile of the membrane and thick dashed lines represent the bound and free filaments. }
\label{model} \end{figure}

The filaments are modeled as rigid rod-like polymers, composed of monomers of length $d$ \cite{kolomeisky2015, Sadhu2016, Sadhu2018, Sadhu2019}. In Fig. \ref{model}, we have represented the filaments using dashed lines, each dash denoting a monomer. A (de)polymerization event increases (decreases) the length of a filament by an amount $d$. For the sake of simplicity, we have used $d =\delta$ here [also see Eq. (\ref{eq:energy_cost})]. In Sec. \ref{sup:sec7} of \cite{sup}, we have included our data for the $d\neq\delta$ case. There are two types of filaments in our model: a free filament, which is not in contact with the membrane, and a bound filament, whose tip is in contact with the membrane site \cite{kirone,Sadhu2016,Sadhu2018,Sadhu2019}. The point of contact is referred to as a binding site. For a free filament, polymerization happens with rate $U_0$. For a  bound filament, however, a polymerization process increases the height of the binding site by an amount $d$ and hence there is an energy cost involved in this process. For a positive (negative) energy cost, the bound filament polymerization rate is taken to be $U_0 R_+$  ($U_0 R_-$), while for zero energy cost, the rate is simply $U_0$. The depolymerization rate of the filament is equal to $W_0$ always, as it does not involve any membrane movement. Note that, in the absence of any polymerization force from the filaments, the membrane tends to stay flat and this aspect is somewhat similar to lamellipodial protrusions, rather than filopodial protrusions observed in a cell. In our simulation, each Monte Carlo step consists of $L$ membrane updates and $N$ filament updates. The detailed simulation algorithm has been presented in Sec. \ref{sup:sec1} of \cite{sup}.

\section{Results}
\label{sec:results}
We present our results in this section. As mentioned in the Introduction, our main result is the non-monotonic variation of the membrane velocity as a function of $\sigma$. We first present our simulation data showing this effect and then we show how this non-trivial behavior can be explained from detailed measurement of the energy cost and the shape of the membrane. We show our results for a single filament ($N=1$) here. Most of our conclusions remain valid even for the case of multiple filaments, if the filament density is not too high.

\subsection{Membrane velocity shows a dip and a peak with $\sigma$}
\label{sec:sigma_v}
Pushed by the growing filaments, the membrane develops an average velocity. In Fig. \ref{sigma_vs_v}, we present our data for the variation of membrane velocity $V$  as a function of the surface tension $\sigma$, for a fixed value of bending rigidity $\kappa$. We find that for small $\kappa$, velocity decreases monotonically with $\sigma$, as expected \cite{Sadhu2018}. However, as $\kappa$ is held fixed at a moderate or large value, $V$ shows a rich behavior: starting with a non-zero value at $\sigma=0$, $V$ first decreases with $\sigma$ and reaches a minimum and then it increases to reach a maximum before finally decreasing exponentially for large $\sigma$. While this non-monotonic behavior is in general interesting \cite{Sadhu2018}, the most intriguing observation here is that, there is a range of $\sigma$ for which $V$ grows with $\sigma$. This growth is counter-intuitive because one generally expects that with increasing surface tension, it should become more difficult for the filament to push the membrane. As we show in the following section, this expectation in fact breaks down.
\begin{figure}[h!]
\centering
\includegraphics[scale=0.7]{sigma_vs_v_diff_kappa_N_1_L_64.eps}
\caption{$\sigma-V$ curve for different values of $\kappa$.  For small $\kappa$ values, $V$ decreases with $\sigma$, but for moderate or large $\kappa$ values, $V$ shows a minimum and a maximum. We scale each curve by $V_{max}$, which denotes the largest value of $V$ for a particular $V-\sigma$ curve. We have $V_{max} = \num{3.45e-2} nm/sec$, $\num{1.55e-3} nm/sec$ and $\num{2.00e-4} nm/sec$ for $\kappa=0.1 pN/nm$, $0.8 pN/nm$ and $1.2 pN/nm$ respectively. The vertical lines show $\sigma$ values used in Fig. \ref{avg_height}. Left inset: peak position $\sigma^*$ shifts leftward  with $\kappa$. Right inset: Exponential fall of $V$ for large $\sigma$ with $\kappa=1.2 pN/nm$. The decay constant $3.6$ matches closely with the analytical prediction $2\beta \delta ^2$. Here, for all the plots, we use $L=64$. The filament depolymerization rate, $W_0 = 1.4 s^{-1}$ \cite{review1, pollard,kirone}, the free filament polymerization rate, $U_0=2.784 s^{-1}$ \cite{review1,pollard,kirone}, the monomer size is $d=\delta=2.7 nm$ \cite{review1, hansda2014,kirone}, $\beta=1/k_B T$, and we have used $T=300 K$.}
\label{sigma_vs_v} \end{figure}

The data in Fig. \ref{sigma_vs_v} is for a fixed $L$ value. We have checked that as $L$ is increased, $V \sim 1/L$, but the peak position does not depend on $L$ for sufficiently large $L$ values. As $L$ becomes small, the peak in the $V-\sigma$ plot shifts to lower $\sigma$ values with decreasing $L$ (Fig. \ref{sigma_vs_v_diff_L} of \cite{sup}).

\subsection{Bending energy cost shows a peak with $\sigma$ }
\label{sec:4th_derivative}
From Eq. (\ref{eq:energy_cost}), it follows that the energy cost $\mathcal{E}_b^+$ for creating a protrusion at the binding site $i=b$ can be decomposed into two parts, the surface energy cost $\Sigma_b^+$ and the bending energy cost $\mathcal{K}_b^+$. In Fig. \ref{sigma_vs_del_e}, we plot the steady state average values of $\mathcal{E}_b^+$, $\Sigma_b^+$ and $\mathcal{K}_b^+$ as a function of $\sigma$, for a large value of $\kappa$. While $\Sigma_b^+$ increases with $\sigma$ monotonically, $\mathcal{K}_b^+$ is found to show a peak. Moreover, this peak appears for relatively small $\sigma$ values, when the total energy cost is actually controlled by the bending interaction. As a result, $\mathcal{E}_b^+$ also shows a peak and immediately after the peak, there is a range of $\sigma$ values, for which energy cost decreases with $\sigma$, contrary to the normal intuition. Finally, for very large $\sigma$, when the membrane is almost flat, $\mathcal{K}_b^+$ saturates and $\Sigma_b^+$ increases linearly with $\sigma$ [see Eq. (\ref{eq:energy_cost})] and $\mathcal{E}_b^+$ also shows a linear rise.
\begin{figure}[h!]
\centering
\includegraphics[scale=0.7]{sigma_vs_del_E_b_k_1_pt_2_N_1_L_64.eps}
\caption{Average energy cost for creating a protrusion at the binding site. Surface energy cost $\Sigma_b^+$ increases with $\sigma$ but bending energy cost $\mathcal{K}_b^+$ shows a peak. For moderate or large $\kappa$ values, when $\mathcal{K}_b^+$ is dominant, this gives rise to a peak in total energy cost $\mathcal{E}_b^+$. For very large $\sigma$, we see  $\mathcal{K}_b^+$ saturating and variation in $\mathcal{E}_b^+$ is then controlled by $\Sigma_b^+$ again. These data are for $\kappa=1.2 pN/nm$ and the other simulation parameters are as in Fig. \ref{sigma_vs_v}. The vertical lines correspond to $\sigma$ values used in Fig. \ref{avg_height}.}
\label{sigma_vs_del_e} \end{figure}

Although, it seems quite surprising that energy cost shows non-monotonic variation with $\sigma$, this effect can be very simply explained from the consideration of the membrane shape near the binding site. For very large $\sigma$ and $\kappa$, the membrane remains almost flat and higher order derivatives of height can be ignored and one can assume an almost linear height variation around the binding site. We check this explicitly in our simulations [see Fig. \ref{avg_height}(a)]. Note that, here we only present the height profile near the binding site and far away from this point, height variation does not remain linear anymore, the profile gradually flattens out. However, this part of the membrane does not influence the energy cost near the binding site. In the region around the binding site, where the height profile is linear, the bending energy, which depends on the second derivative of height, has a non-zero value only at the binding site, and zero elsewhere. Keeping $\kappa$ fixed at large value, as we lower $\sigma$,  the qualitative shape of the membrane remains similar, but the magnitude of the height gradient is now larger [Fig. \ref{avg_height}(b)]. This, however, gives rise to a very high bending energy around the binding site, since height gradient changes sharply from a large positive value to a large negative value across the binding site. Therefore, as $\sigma$ is lowered further, such a membrane shape becomes unsustainable and instead a shape shown in Fig. \ref{avg_height}(c) is observed, where the height gradient does not change so sharply but varies more slowly in space. Such a membrane shape stores less bending energy than the one with a linear height variation. In Sec. \ref{sup:sec3} of \cite{sup}, we have also included an approximate analytical explanation of this effect, where we analytically calculated the membrane shape, adapting a simplified description outlined in \cite{nirgov08}.
\begin{figure}[h!]
\centering
\includegraphics[scale=0.7]{avg_height_4th_deri_large_k.eps}
\caption{(a)-(c): Local shape of the membrane around the binding site for different $\sigma$ with $\kappa=1.2 pN/nm$. (d): Discrete fourth derivative of height at the binding site, $h_b''''=h_{b-2}-4h_{b-1}+h_b-4h_{b+1}+h_{b+2}$, as a function of $\sigma$. Here, $L=64$ and the other simulation parameters are as in Fig. \ref{sigma_vs_v}.}
\label{avg_height} \end{figure}

With the above picture in mind, it is now easy to understand the variation of the fourth derivative of height (that determines the bending energy required to create a protrusion) at the binding site. In the configuration shown in Fig. \ref{avg_height}(a), the second derivative has a negative value at the binding site and zero elsewhere in its neighborhood. As we lower $\sigma$, the configuration changes to what is shown in Fig. \ref{avg_height}(b) and the magnitude of the second derivative at the binding site  increases, {\sl i.e.} the minimum in $h_i^{\prime \prime}$ is now sharper. This will increase its curvature, {\sl i.e.} the fourth derivative at the binding site. However, as the membrane shape changes to Fig. \ref{avg_height}(c) for even lower $\sigma$, $h_i^{\prime \prime}$ now has a non-zero (negative) value even some distance away from the binding site, since the height profile is not linear anymore. Around the binding site, $h_i^{\prime \prime}$ now varies more gradually and while it still has a minimum at the binding site, the minimum is not as sharp. This corresponds to a lower curvature and the fourth derivative decreases in magnitude. We show this explicitly in Fig. \ref{avg_height}(d).

However, when $\kappa$ is small, the above effect is absent. Since the bending energy is actually responsible for the change in membrane shape and for small $\kappa$, the bending energy is just not significant enough to bring about this change. In this case, the membrane height profile around the binding site remains linear for all $\sigma$ values and the fourth derivative falls monotonically with $\sigma$ \cite{sup}. As a result, energy cost and velocity varies monotonically with $\sigma$ (see Figs. \ref{shape_small_k} and \ref{kappa_vs_v} in \cite{sup}).

\subsection{A more quantitative explanation of the full $V-\sigma$ curve}
\label{sec:exp}

In the previous section, we have offered a simple qualitative explanation of how the energy cost for creating a protrusion shows a peak with $\sigma$, for large values of $\kappa$. This makes it plausible that $V$ can also show a similar peak, since pushing the membrane upward becomes easier as $\sigma$ is increased within a certain range. However, the net membrane velocity results from both upward and downward movements of the membrane and in this section, we present a more quantitative and detailed calculation, that explains the full $V-\sigma$ curve, along with its maximum and minimum.

The membrane being in steady state, its average velocity should be the same at all sites and below, we write down the expression for velocity at the binding site $b$. For large $\kappa$, height gradients at the binding site are small and
hence energy costs $\mathcal{E}_b^ \pm$ are positive [also see Eq. (\ref{eq:energy_cost})], which gives
\begin{equation}
V_b=\delta [U_0 p_0 e^{-\beta \mathcal{E}_b^+} + e^{-\beta \mathcal{E}_b^+} - (1-p_0) e^{-\beta \mathcal{E}_b^-}],
\label{eq:vb} \end{equation}
where $p_0$ is the contact probability of the filament tip with the  membrane. The first term in Eq. (\ref{eq:vb}) corresponds to the bound filament polymerization. The second and third terms represent thermal fluctuation of the membrane height at the binding site. Here, we have used the fact that $h_b$ can always increase with rate $e^{-\beta {\mathcal E}_b^+}$, but it can decrease only when the filament is not bound to the membrane.

As explained in Eq. (\ref{eq:energy_cost}), $\mathcal{E}_b^+$ ($\mathcal{E}_b^-$) is the energy cost for increasing (decreasing) $h_b$ by an amount $\delta$, and generally it depends on the local configuration around the binding site. Here, we use the approximation that $\mathcal{E}_b^\pm$ may be replaced by their average values (to keep our notations simpler, we have used the same symbol for the average quantities as well). We have also checked that (data not shown here) as long as $\kappa$ is not too small, the contact probability $p_0$ does not show much variation and remains close to $1/2$, its value for a rigid barrier \cite{Sadhu2018}.

For large $\kappa$ and small $\sigma$, the surface energy cost is much smaller than the bending energy cost and we can replace $\mathcal{E}_b^\pm \approx \mathcal{K}_b^\pm$ and this gives
\begin{equation}
V \approx e^{-6 \beta \delta ^2 \kappa} \delta \{(U_0+1)p_0 e ^{-2 \beta \delta \kappa h_b''''}-2(1-p_0) Sinh{(2\beta \delta \kappa h_b'''')}\}.
\label{eq:vk} \end{equation}

It is easy to see from the above expression that as $h_b^{\prime \prime \prime \prime}$ first increases with $\sigma$, reaches a peak and then decreases, $V$ also initially decreases with $\sigma$, reaches a minimum, and then starts increasing. As $\sigma$ increases further, it is not possible to neglect the elastic energy anymore and Eq. (\ref{eq:vk}) does not remain valid. However, in this case, both $\sigma$ and $\kappa$ are large, and the magnitudes of $h_b^{\prime \prime}$ and $h_b^{\prime \prime \prime \prime}$ are negligible, and we can write
\begin{equation}
V \simeq \delta  (U_0+1) p_0 e^{-6\beta \delta ^2 \kappa} e^{-2 \beta \delta ^2  \sigma}, \label{eq:vks}
\end{equation}
where $V$ decreases exponentially with $\sigma$. We verify this from our numerical simulation (Fig. \ref{sigma_vs_v}, right inset). The exponential decay constant numerically observed is $3.6$, which is close to the analytically predicted value $2\beta \delta ^2 \simeq 3.52$. When $\kappa$ is held fixed at a larger value, the exponential decay in Eq. (\ref{eq:vks}) starts at a smaller value of $\sigma$, since the height derivatives become negligible already. This is why the peak position $\sigma^\ast$ shifts towards smaller values as $\kappa$ increases (Fig. \ref{sigma_vs_v}, left inset).

\section{Discussions}
\label{sec:sum}
Throughout this rapid communication, we have limited our studies to one dimension and it is an important question whether our conclusions remain valid for a two dimensional membrane as well. Generalizing Eqs. (\ref{eq:hamiltonian}) and (\ref{eq:energy_cost}) for the two dimensional case, it can be shown that the bending energy cost is much higher in this case, which slows down the time evolution for moderate or high values of $\kappa$. Although, for small $\kappa$ values we have been able to verify (data not shown here) that $V$ decreases monotonically with $\sigma$, as found in the one dimensional system, higher $\kappa$ values remain numerically inaccessible to us. More research is needed to conclude with certainty if non-monotonic variation of $V$ with $\sigma$ for high $\kappa$ values persists in the two dimensional model as well.

It should be possible to experimentally verify our conclusions. One direct and perhaps the simplest measurement would be to measure the membrane velocity for different values of $\sigma$ and $\kappa$ \cite{marcy2004, baudry2011, theriot2005} and see whether a non-monotonicity as shown in Fig. \ref{sigma_vs_v} can be found. The energy scales involved in our simulations are actually comparable to physical systems. For example, the increase of $V$ with $\sigma$ that we observe for high $\kappa$, corresponds to the surface energy varying in the range of $0.5 $  to $2$ ($ 10^{-19}$ $J$) and the bending energy lies in between $1$ and $3.5$ ($ 10^{-19}$ $J$). These values are within the experimentally observed ranges for flexible membranes \cite{Dimova}. Although, in real systems, many other factors, apart from surface energy-bending energy interplay, are relevant, it would be interesting to see whether this basic signature of a protrusion formation mechanism can still be found.

\section{Acknowledgements}
We acknowledge useful discussions with P. Pradhan and A. Kundu. SC acknowledges financial support from the Science and Engineering Research Board, India (Grant No. EMR/2016/001663). The computational facility used in this work was provided through the Thematic Unit of Excellence on Computational Materials Science, funded by Nanomission, Department of Science and technology (India).


\section{Introduction}
\label{sec:intro}
Image restoration is the task of recovering a sharp image $u$ starting from a blurred and noisy observation $g$. In this work, we consider a degradation model of the form
\begin{equation}
\label{eq:degmod}
g = Ku + e\,,
\end{equation}
where $g,u\in\mathbb{R}^n$ are vectorized images, $K \in \mathbb{R}^{n\times n}$ is the linear blur operator and $e\in\mathbb{R}^n$ is an additive noise vector. A possible strategy to overcome the ill-posedeness of the linear system in \eqref{eq:degmod} is to reformulate the problem. We thus look for $u^*$, estimate of the original $u$, which solves a well-posed problem. In the \emph{variational} approach, $u^*$ is the minimizer of a cost functional $\mathcal{J}(u;K,g):\mathbb{R}^n\to \mathbb{R}$. In formula,
\begin{equation}
\label{eq:varmod}
u^*\in\arg\min_{u\in\mathbb{R}^n} \; \big\{\mathcal{J}(u;K,g):=\mathcal{R}(u) + \mu \mathcal{F}(u;K,g)\big\}.
\end{equation}
The functionals $\mathcal{R}$ and $\mathcal{F}$ are commonly referred to as the \emph{regularization} and the \emph{data fidelity} term, respectively. While $\mathcal{R}$ encodes prior information on the desired image $u$, $\mathcal{F}$ is a data term which measures the `distance' between the given image $g$ and $u$ after the action of the operator $K$ with respect to some norm corresponding to the noise statistics in the data, cf., e.g., \cite{stuart}. The regularization parameter $\mu > 0$ controls the trade-off between the two terms.
\\In this paper, we consider an Additive White Gaussian Noise (AWGN) corrupting the blurred image $Ku$, i.e. $e\sim\mathcal{N}(0,\sigma^2I_n)$, where $I_n$ is the $n$-dimensional identity matrix.
It is well known that, in presence of AWGN, a suitable choice for $\mathcal{F}(u;K,g)$ is the so-called L$_2$ fidelity term, reading as,
\begin{equation}
\label{eq:l2}
\mathcal{F}(u;K,g) = \mathrm{L}_2(u;K,g) = \frac{1}{2} \lVert Ku - g \rVert_2^2.
\end{equation}
A popular choice for the regularization term $\mathcal{R}(u)$ is given by the TV semi-norm \cite{ROF},
\begin{equation}
\mathcal{R}(u) = \mathrm{TV}(u)=\sum_{i=1}^{n}  \lVert (\nabla u)_{i} \rVert_2 \, ,
\label{eq:TV}
\end{equation}
where $(\nabla u)_i := \big( (D_h u)_i , (D_v u)_i \big)^T \in \mathbb{R}^2$ denotes the discrete gradient of image $u$ at pixel $i$,
with $D_h,D_v \in \mathbb{R}^{n \times n}$ linear operators representing finite difference discretizations of the first-order horizontal and vertical partial derivatives, respectively.
\\Coupling the L$_2$ data term with the TV regularizer leads to one of the most widespread variational models for image restoration problem, the TV-L$_2$ (or ROF) model,
\begin{equation}
\label{eq:tvrof}
u^*\in\arg\min_{u\in\mathbb{R}^n} \bigg\{\sum_{i=1}^{n}  \lVert (\nabla u)_{i} \rVert_2 +
\frac{\mu}{2} \lVert Ku - g \rVert_2^2\bigg\}.
\end{equation}
The global perspective of the TV-L$_2$ model does not allow to diversify the action of the regularizer on regions of the image presenting different properties. In \cite{siam,ncmip,cmbbe,vip}, the authors have proposed space-variant regularization term based on statistical assumptions on the distribution of the $\ell_2$-norm of the gradients and on the gradients themselves.
\\In this paper, we propose a space-variant anisotropic extension of the TV regularizer in (\ref{eq:TV}) which, as it will be illustrated in Sect. \ref{sec:map}, comes from the \emph{a priori} assumption that the gradients of the target image $u$ distribute locally according to a Bivariate Laplace Distribution (BLD). The proposed BLTV regularizer takes the form
\begin{eqnarray}
\label{eq:jlapreg}
\mathrm{BLTV}(u;\lambda_1,\lambda_2,\theta) &=& \sum_{i=1}^{n} \lVert \Lambda_i R_{\theta_i} (\nabla u)_i\rVert_1\\
&=& \sum_{i=1}^{n}\bigg[\lambda_i^{(1)}|\langle r_i^{(1)},(\nabla u)_i \rangle| + \lambda_i^{(2)} |\langle r_i^{(2)},(\nabla u)_i \rangle|\bigg]\,,
\end{eqnarray}
where $\Lambda_i$ is a $2\times 2$ positive definite diagonal matrix and $R_{\theta_i}$ is the rotation matrix corresponding to the angle $-\theta_i$. Mathematically,
\begin{equation}
\label{eq:lamrth}
\Lambda_i = \begin{pmatrix}
\lambda_i^{(1)} & 0 \\
0 &\lambda_i^{(2)}
\end{pmatrix}\,,\quad
R_{\theta_i} = \begin{pmatrix}
\cos \theta_i & \sin \theta_i \\
-\sin \theta_i &\cos \theta_i
\end{pmatrix} =
\begin{pmatrix}
r_i^{(1)}\\
r_i^{(2)}
\end{pmatrix}.
\end{equation}
We denote by $\lambda_1,\lambda_2,\theta \in \mathbb{R}^n$ the maps of the parameters defining the local distributions. Hence, the proposed BLTV-L$_2$ variational restoration model reads as
\begin{equation}
\label{eq:jtvl2}
u^*\in\arg\min_{u\in\mathbb{R}^n} \bigg\{\mathrm{BLTV}(u;\lambda_1,\lambda_2,\theta)+
\frac{\mu}{2} \lVert Ku - g \rVert_2^2\bigg\}.
\end{equation}
The $3n$ free parameters defining the BLTV regularizer in (\ref{eq:jlapreg})--(\ref{eq:lamrth}) hold the potential for effectively modelling local image properties, in particular driving in a suitable adaptive manner the strength and direction of non-linear TV-diffusion.
In Figs.\ref{fig:pdf}(a)-(b) the red ellipses represent  TV-diffusion strengths along all possible directions at few sample pixel locations for TV and BLTV regularizers. It is evident how for classical TV such ellipses turn out to be circles (isotropy) of constant radius (space-invariance), whereas our BLTV regularizer allows for  ellipses (anistropy) of different size (space-variance).
In practice, such flexibility of BLTV can be (and will be) exploited to diffuse in different ways in regions exhibiting different properties: for instance, strong isotropic diffusion in homogeneous regions, strongly anisotropic diffusion in regions characterized by a dominant edge direction.
In Figs.\ref{fig:pdf} (c)-(d) we show the level curves of the TV regularizer and one among the infinity of possible configuration of the level curves of the BLTV regularizer, respectively, revealing once again the flexibility of the proposed regularizer.

\begin{figure}[tbh]
	\center
	\begin{tabular}{cccc}
		\includegraphics[width=1.1in]{images/diffTV.pdf} &
		\includegraphics[width=1.1in]{images/diffBTV.pdf}&
			\includegraphics[width=1.1in]{images/curveTV.pdf} &
		\includegraphics[width=1.1in]{images/curveBTV.pdf}\\
		(a)&(b)&(c)&(d)\\
	\end{tabular}
	\caption{Space-invariant isotropic TV-diffusion (a), space-variant anisotropic BLTV-diffusion (b), level curves of TV (c) and BLTV (d) regularization terms.}
	\label{fig:pdf}
\end{figure}

Together with the variational model in \eqref{eq:jtvl2}, we also propose an efficient and robust procedure for automatically estimating the parameter maps $\lambda^{(1)}$, $\lambda^{(2)}$, $\theta$ based on a Maximum Likelihood (ML) approach. The parameter maps can be updated along the iterations of the Alternating Direction Methods of Multipliers (ADMM), which is the algorithm adopted here to solve the minimization problem. Notice that the convexity of the BLTV-L$_2$ model ensures the convergence of the ADMM.

\vspace{-0.5cm}

\section{Deriving the model via MAP}
\label{sec:map}
Applying the \emph{Maximum A Posteriori} (MAP) estimation approach to image restoration consists in computing the restored image as a global maximizer of the posterior probability Pr$(u|g;K)$ of the unknown target image $u$ given the observation $g$, namely:
\begin{equation}
u^* \:\;{\in}\;\: \arg \max_{u \in \mathbb{R}^n}\;
\text{Pr}(u|g;K) = \arg \min_{u \in \mathbb{R}^n} \;
\left\{ \,-\mathrm{ln} \text{Pr}(g|u;K)-\mathrm{ln} \text{Pr}(u)
\, \right\}\,,
\label{eq:map1}
\end{equation}
where, after applying the Bayes' rule, we drop the \emph{evidence} term Pr$(g)$ and extracted $-\ln$ of the objective function. The two terms Pr$(u)$ and Pr$(g|u;K)$ in
(\ref{eq:map1}) are referred to as the \emph{prior} and the \emph{likelihood}.
The likelihood term associated with AWGN corruption takes the form
\begin{equation}
\mathrm{Pr}(g|u;K)=\prod_{i=1}^{n}\,\frac{1}{\sqrt{2\pi}\sigma}\,\text{exp}\bigg(-\frac{(Ku-g)_{i}^{2}}{2\sigma^{2}}\,\bigg)=W\,\text{exp}\bigg(-\frac{\lVert Ku-g\rVert_{2}^{2}}{2\sigma^{2}}\,\bigg),
\label{eq:g_like}
\end{equation}
where $\sigma>0$ denotes the AWGN standard deviation and $W > 0$ is a normalization constant.
For what concerns the prior, a common choice is to model the unknown image $u$ as a Markov Random Field (MRF) such that the image can be characterized by its Gibbs prior distribution, whose general form is:
\begin{equation}
\text{Pr}(u)= \prod_{i = 1}^{n} \text{Pr}_i(u) = \prod_{i = 1}^{n} z_i \, \text{exp}\,(\,-\alpha \, V_{c_{i}}(u)\,)=  Z \, \text{exp}\,\bigg(\,-\alpha\,\sum_{i = 1}^{n}  V_{c_{i}}(u)\,\bigg),
\label{eq:mrf}
\end{equation}
where $\alpha>0$ is the MRF parameter, $\{c_{i}\}_{i=1}^{n}$ is the set of all cliques (a clique is a set
of neighboring pixels) for the MRF, $V_{c_{i}}$ is the potential function defined on the clique $c_{i}$
and $Z = \prod_{i = 1}^{n} z_i$ is a normalization constant. Choosing as potential function at the generic $i$-th pixel the magnitude of the discrete gradient at the same pixel, i.e. $V_{c_{i}}=\lVert (\nabla u)_{i} \rVert_{2}$ for any $i \in \{1,\ldots,n\}$, the Gibbs prior in \eqref{eq:mrf} reduces to the TV prior which, plugged into (\ref{eq:map1}), yields the popular TV regularizer. We remark that the TV prior corresponds to choosing
\begin{equation}
\label{eq:Pri_TV}
\text{Pr}_i(u) = z_i \, \text{exp}\,\left(\,-\alpha \lVert (\nabla u)_{i} \rVert_{2} \right)\, .
\end{equation}
Here, we propose a space-variant anisotropic generalization of the local prior $\text{Pr}_i(u)$ above. More in detail, we assume that the discrete gradient at any pixel $i$ distributes according to a space-variant BLD, such that our prior reads as
\begin{eqnarray}
\notag
\mathrm{Pr}(u) &=& \prod_{i = 1}^{n} \text{Pr}_i(u) =
\prod_{i = 1}^{n} z_i \,
\exp\bigg(- \lVert \Lambda_iR_{\theta_i} (\nabla u)_i \rVert_1\bigg)\\ &=& Z \, \exp\bigg(-\sum_{i=1}^{n}\lVert \Lambda_iR_{\theta_i} (\nabla u)_i \rVert_1\bigg)\,,
\label{eq:jtvpr}
\end{eqnarray}
where $Z$ is the normalization constant and $\Lambda_i, R_{{\theta}_i} \in \mathbb{R}^{2 \times 2}$ are defined in (\ref{eq:lamrth}).

Plugging the prior \eqref{eq:jtvpr} and the likelihood \eqref{eq:g_like} into \eqref{eq:map1} and neglecting the  constant terms, the proposed BLTV-L$_2$ variational model \eqref{eq:jtvl2} is obtained.

\vspace{-0.5cm}

\section{Parameter estimation via ML}
\label{sec:ml}
In order to make the introduction of the proposed regularizer actually useful, an efficient, robust, and automatic procedure for the estimation of the  parameter maps $\lambda^{(1)},\lambda^{(2)},\theta\in \mathbb{R}^n$ identifying all the local BLDs has to be proposed as well. To this aim, we resort to the ML approach. Consider a set of $N$ 2-dimensional samples
$\mathcal{S} := \left\{ s_1,\ldots,s_N \right\}$ drawn from a BLD with parameters $(\lambda_1,\lambda_2,\theta)$. Here, the samples play the role of image gradients at pixels of a neighborhood of radius $r$ centered at a generic pixel $i$.
Assuming independence of the samples, according to the definition of BLD, the likelihood function is defined by
\begin{equation}
 \label{eq:ml}
 \mathrm{Pr}(\mathcal{S}|\lambda_1,\lambda_2,\theta) = \bigg(\frac{\lambda_1\,\lambda_2}{4}\bigg)^N \exp\bigg(-\sum_{i=1}^{N}\bigg(\lambda_1|\langle r_1, s_i\rangle|+\lambda_2|\langle r_2, s_i\rangle|\bigg)\bigg)\,,
 \end{equation}
where, clearly, $r_1$ and $r_2$ depend on $\theta$ - see \eqref{eq:lamrth}.
The goal here is to find
\begin{eqnarray}
\notag
(\lambda_1^*,\lambda_2^*,\theta^*) &\in& \arg\max_{\lambda_1,\lambda_2,\theta} \mathrm{Pr}(\mathcal{S}|\lambda_1,\lambda_2,\theta) = \arg\min_{\lambda_1,\lambda_2,\theta}  -\mathrm{ln} \mathrm{Pr}(\mathcal{S}|\lambda_1,\lambda_2,\theta)\\
&=&-N\mathrm{ln}\frac{\lambda_1\,\lambda_2}{4} + \sum_{i = 1}^{N}\bigg(\lambda_1|\langle r_1, s_i\rangle|+\lambda_2|\langle r_2, s_i\rangle|\bigg).
\label{eq:obj_ml}
\end{eqnarray}
Imposing a first-order optimality condition with respect to $\lambda_1,\lambda_2$ leads to the following closed-form estimation formulas:
\begin{equation}
\label{eq:updm}
\lambda_1 = \bigg(\frac{1}{N}\sum_{i = 1}^{N}|\langle r_1, s_i\rangle|\bigg)^{-1}\,,\quad
\lambda_2 = \bigg(\frac{1}{N}\sum_{i = 1}^{N}|\langle r_2, s_i\rangle|\bigg)^{-1}.
\end{equation}
Substituting the expressions in \eqref{eq:updm} into the objective function \eqref{eq:obj_ml}, we thus obtain the simplified minimization problem in the only variable $\theta$:
\begin{equation}  \label{pb:min_theta}
\theta^*\in \arg \min_{\theta}\bigg\{\mathrm{ln}\bigg(\sum_{i=1}^{N}|\langle r_1, s_i\rangle|\bigg)+\mathrm{ln}\bigg(\sum_{i=1}^{N}|\langle r_2, s_i\rangle|\bigg)\bigg\}\,.
\end{equation}

\vspace{-0.5cm}

\section{ADMM}
\label{sec:admm}
In order to solve numerically the proposed image restoration model \eqref{eq:jtvl2}, we use an ADMM-based algorithm - see \cite{boyd}. We first introduce two auxiliary variables $w \in \mathbb{R}^n$, $t \in \mathbb{R}^{2n}$ and rewrite the model in the equivalent linearly constrained form:

\vspace{-0.65cm}
\begin{eqnarray}\label{eq:PM_ADMM_a}
\{ \, u^*,w^*,t^* \}&
\:\;{\in}\;\:&
\arg\min_{u,w,t}
\bigg\{ \:
\sum_{i = 1}^{n}\lVert \Lambda_i R_{\theta_i} t_i\rVert_1
\;{+}\;
\frac{\mu}{2} \, \| w \|_2^2
\: \bigg\}\\
&&\mathrm{subject}\;\mathrm{to}:\quad
w \;{=}\; K u - g, \;\:
t \;{=}\; Du ,
\end{eqnarray}
where the space-variant matrices $\Lambda_i$, $R_{\theta_i}$ can be estimated via the ML procedure described in Sect. \ref{sec:ml} based only on the observed image $g$ (i.e. as a preliminary pre-processing step) or also updated along the ADMM iterations.
We define the augmented Lagrangian functional:
\begin{eqnarray}\nonumber
\mathcal{L}(u,w,t;\rho_w,\rho_t)&:=& \sum_{i = 1}^{n}\lVert \Lambda_i R_{\theta_i} t_i \rVert_1+\frac{\mu}{2} \| w \|_2^2 - \rho_t^{T}(t - D u) +\frac{\beta_t}{2} \| t - D u \|_2^2\\
\label{eq:PM_AL}
&-&\rho_w^{T}(w - (Ku-g))+ \frac{\beta_w}{2}\| w - (Ku-g) \|_2^2,
\end{eqnarray}
where $\beta_w, \beta_t > 0$ are scalar penalty parameters and $\rho_w \in \mathbb{R}^n$, $\rho_t \in \mathbb{R}^{2n}$
are the vectors of Lagrange multipliers associated with the given linear constraints.
The solution $\{u^*,w^*,t^*\}$ of problem (\ref{eq:PM_ADMM_a}) is a saddle point for $\mathcal{L}$ in (\ref{eq:PM_AL}),
{see, e.g., \cite{boyd}}. Hence, we can alternate a minimization step with respect to $t,u,w$ with a maximization step with respect to $\rho_t,\rho_w$. Mathematically,

\begin{eqnarray}
&
u^{(k+1)} &
\;{\leftarrow}\;\;\,\,
\arg\min_{u \in \mathbb{R}^n} \;
\mathcal{L}(u,w^{(k)},t^{(k)};\rho_w^{(k)},\rho_t^{(k)}) \, ,
\label{eq:PM_ADMM_u} \\
&
w^{(k+1)} &
\;{\leftarrow}\;\;\,\,
\arg\min_{r \in \mathbb{R}^n} \;
\mathcal{L}(u^{(k+1)},w,t^{(k)};\rho_w^{(k)},\rho_t^{(k)}) \, ,
\label{eq:PM_ADMM_r} \\
&
t^{(k+1)} &
\;{\leftarrow}\;\;\,\,
\arg\min_{t \in \mathbb{R}^{2n}} \;
\mathcal{L}(u^{(k+1)},w^{(k+1)},t;\rho_w^{(k)},\rho_t^{(k)}) \, ,
\label{eq:PM_ADMM_t} \\
&
\rho_w^{(k+1)} &
\;{\leftarrow}\;\;\,\,
\rho_w^{(k)} \;{-}\; \beta_r \, \big( \, w^{(k+1)} \;{-}\; (K u^{(k+1)}-g) \, \big) \, ,
\label{eq:PM_ADMM_lr} \\
&
\rho_t^{(k+1)} &
\;{\leftarrow}\;\;\,\,
\rho_t^{(k)} \;{-}\; \beta_t \, \big( \, t^{(k+1)} \;{-}\; D u^{(k+1)} \, \big) \, .
\label{eq:PM_ADMM_lt}
\end{eqnarray}

\noindent {The solution of the primal sub-problem \eqref{eq:PM_ADMM_u} can be efficiently computed by means of standard linear Fast Fourier Transform (FFT) solvers. The sub-problem \eqref{eq:PM_ADMM_r} can be solved in closed-form by following \cite[Section 3]{ncmip}. Finally, the sub-problem \eqref{eq:PM_ADMM_t} can be solved by computing efficiently the proximal operator of the anisotropic $1$ norm, for which the proof of \cite[Proposition 6.3]{siam} can be easily adapted. The regularization parameter $\mu$ is updated along the iterations so as to fulfill the
global discrepancy principle as described in \cite{ape}. We refer the reader also to \cite{cmbbe,tvp} for more details on the numerical solution of the algorithm.}

\vspace{-0.5cm}

\section{Experimental results}
\label{sec:test}
In this section, we evaluate the performance of the proposed BLTV-L$_2$ restoration model compared with the baseline TV-L$_2$ model, also solved by ADMM. The stopping criteria of the ADMM for both models are defined based on the number of iterations as well as on the iterates relative change, i.e. we stop iterating as soon as
\begin{equation}
k \geq 1500 \;\;\; \mathrm{or}\;\;\; \delta^{(k)}:= \frac{\lVert u^{(k)} - u^{(k-1)} \rVert}{\lVert u^{(k-1)}\rVert} \geq 10^{-6}.
\end{equation}
The quality of the restored images $u^*$ is measured by means of the Improved Signal-to-Noise Ratio $\mathrm{ISNR}(u^*,g,u) = 10\log_{10}\frac{\|g-u\|_2^2}{\|u^*-u\|_2^2}$, with $u$ denoting the original uncorrupted image, and of the Structural-Similarity-Index (SSIM) \cite{ssim}.\\

We consider the test image \texttt{brain} in Fig.\ref{fig:data} (a) (570 $\times$ 430) and the test image \texttt{abdomen} in Fig.\ref{fig:data2} (a) (350 $\times$ 480) with pixel values between 0 and 255, synthetically corrupted by space-invariant  blur with Gaussian kernel of parameters \verb|band| = 9, \verb|sigma| = 2, and by AWGN of different levels $\sigma \in \{10,20\}$ - see, e.g., Fig.\ref{fig:data} (b) and Fig.\ref{fig:data2} (b). The parameter maps are computed at the beginning starting from the observed image $g$ and then updated every 300 iterations based on the current iterate. The radius of the neighborhoods used for the local parameter estimation has been set equal to 8 and 5 for the \texttt{brain} and \texttt{abdomen} test images, respectively. The reconstructions of \texttt{brain} for $\sigma=20$ and of \texttt{abdomen} for $\sigma=10$ via TV-L$_2$ and BLTV-L$_2$ models are shown in Figs.\ref{fig:data}(c)-(d) and Figs.\ref{fig:data2}(c)-(d), respectively.

\begin{figure}[tbh]
	\center
	\begin{tabular}{cccc}
		\includegraphics[width=1.1in]{images/true_rec.png} &
		\includegraphics[width=1.1in]{images/data_rec.png} &
		\includegraphics[width=1.1in]{images/tv_rec.png} &
		\includegraphics[width=1.1in]{images/jtv_rec.png}\\
		(a)&(b)&(c)&(d)\\
		\includegraphics[width=1.1in]{images/true_zoom.png} &
		\includegraphics[width=1.1in]{images/data_zoom.png} &
		\includegraphics[width=1.1in]{images/tv_zoom.png} &
		\includegraphics[width=1.1in]{images/jtv_zoom.png}\\
		(e)&(f)&(g)&(h)\\
	\end{tabular}
	\caption{\emph{First row}. Original test image \texttt{brain} (a), observed image corrupted by Gaussian blur and AWGN with $\sigma=20$ (b), TV reconstruction (c), BLTV reconstruction (d). \emph{Second row}. Close-up(s) of the first row.}
	\label{fig:data}
\end{figure}
From a visual inspection, the restoration via BLTV-L$_2$ seems to be more neat and less cartooned than the TV reconstructions. As reported in Tab.\ref{tab:1}, the ISNR and SSIM values for the two test images and for different noise levels obtained by the BLTV-L$_2$ model outperform the ones reached by the TV-L$_2$ model.\\
The final parameter maps computed by BLTV-L$_2$ are shown in Figs. \ref{fig:maps}-\ref{fig:maps2}.

\begin{table}[!t]
	\vspace{0.3cm}
	\centering
	\begin{tabular}{l|cc|cc|cc|cc}
	\hline\hline

&\multicolumn{4}{c}{\textbf{\texttt{brain}}}&\multicolumn{4}{c}{\textbf{\texttt{abdomen}}}\\
		\hline\hline
		&\multicolumn{2}{c}{$\sigma=10$}&\multicolumn{2}{c}{$\sigma=20$}&\multicolumn{2}{c}{$\sigma=10$}&\multicolumn{2}{c}{$\sigma=20$}\\
		\hline\noalign{\smallskip}
		&TV-L$_2$&BLTV-L$_2$&TV-L$_2$&BLTV-L$_2$&TV-L$_2$&BLTV-L$_2$&TV-L$_2$&BLTV-L$_2$\\
		\hline
		ISNR&4.27&\textbf{5.52}&5.00&\textbf{6.67}&3.76&\textbf{4.66}&6.10&\textbf{6.77}\\
		SSIM&0.87&\textbf{0.88}&0.83&\textbf{0.85}&0.78&\textbf{0.80}&0.74&\textbf{0.76}\\
		\noalign{\smallskip}\hline\noalign{\smallskip}
	\end{tabular}
	\caption{Maximum ISNR/SSIM values achieved by TV-L$_2$ and BLTV-L$_2$ on \texttt{brain} and \texttt{abdomen} test images corrupted by AWGN of two different levels.}
	\label{tab:1}
\end{table}
\begin{figure}[tbh]
	\center
	\begin{tabular}{ccc}
		\includegraphics[width=1.45in]{images/brain_mu1.png} &
		\includegraphics[width=1.45in]{images/brain_mu2.png} &
		\includegraphics[width=1.45in]{images/brain_theta.png}\\
		$\lambda^{(1)}$&$\lambda^{(2)}$&$\theta$\\
		\end{tabular}
	\caption{Final parameter maps for \texttt{brain} test image.}
	\label{fig:maps}
\end{figure}

\begin{figure}[tbh]
	\center
	\begin{tabular}{cccc}
		\includegraphics[width=1.1in]{images/true_rec_abd.png} &
		\includegraphics[width=1.1in]{images/data_rec_abd.png} &
		\includegraphics[width=1.1in]{images/tv_rec_abd.png} &
		\includegraphics[width=1.1in]{images/jtv_rec_abd.png}\\
		(a)&(b)&(c)&(d)\\
		\includegraphics[width=1.1in]{images/true_zoom_abd.png} &
		\includegraphics[width=1.1in]{images/data_zoom_abd.png} &
		\includegraphics[width=1.1in]{images/tv_zoom_abd.png} &
		\includegraphics[width=1.1in]{images/jtv_zoom_abd.png}\\
		(e)&(f)&(g)&(h)\\
	\end{tabular}
	\caption{\emph{First row}. Original test image \texttt{abdomen} (a), observed image corrupted by Gaussian blur and AWGN with $\sigma=10$ (b), TV reconstruction (c), BLTV reconstruction (d). \emph{Second row}. Close-up(s) of the first row.}
	\label{fig:data2}
\end{figure}

\begin{figure}[tbh]
	\center
	\begin{tabular}{ccc}
		\includegraphics[width=1.45in]{images/abd_mu1.png} &
		\includegraphics[width=1.45in]{images/abd_mu2.png} &
		\includegraphics[width=1.45in]{images/abd_theta.png}\\
		$\lambda^{(1)}$&$\lambda^{(2)}$&$\theta$\\
	\end{tabular}
	\caption{Final parameter maps for \texttt{abdomen} test image.}
	\label{fig:maps2}
\end{figure}

\emph{Computational times.} We tested the joint parameter estimation + reconstruction model on a standard laptop with inbuilt MATLAB software, version 2016b.
As far as the ML parameter estimation of the parameter maps procedure is concerned,
 we notice that the update \eqref{eq:updm} is explicit, thus very cheap, whereas the computation of $\theta^*$ in \eqref{pb:min_theta} requires the solution of an optimisation problem. We solve the problem by line-searching upon a suitable discretization of the parameter space. For Fig. \ref{fig:data}, the ML parameter estimation procedure took 8.45 secs.

The ADMM algorithmic sub-steps with automatic parameter update every 300 iterations computes the numerical solution in 381 secs for the high-resolution image in Fig. \ref{fig:data}. A possible way to accelerate the speed of the algorithm would be the computation the parameter maps only in terms of the given image and not along the iterations, although of course that would render a less accurate result.

\color{black}

\vspace{-0.5cm}

\section{Conclusions and future works}

We presented a new space-variant anisotropic regularization term for image restoration based on the \emph{a priori} statistical assumption that the gradients of the unknown target image distribute locally according to space-variant bivariate Laplace distributions. The high flexibility of the proposed regularizer together with the presented ML parameters estimation procedure and ADMM-based minimization algorithm yield a very effective and efficient approach. Preliminary experiments on images corrupted by blur and AWGN strongly indicate that the proposed variational model achieves high-quality restorations and, in particular, outperforms by far the results obtained by classical TV-regularized restoration models.
Coupling the proposed regularizer with other fidelity terms suitable for dealing with noises other than Gaussian - such as, e.g., Laplace, Poisson and mixed Poisson-Gaussian \cite{MIXED} - is a matter being studied.

\vspace{-0.5cm}

\input{ref}
\end{document}


\section{Introduction}

High power fiber amplifiers have become a widely used tool in research and industry \cite{jauregui} due to their robust operation, excellent beam quality and power scalability. Some applications directly pose high requirements on the stability and the intensity noise characteristics of fiber amplifiers, such as their potential usage as laser sources in next-generation gravitational wave detectors  \cite{kracht}. Other applications, such as laser cutting, require low noise more indirectly, since single mode power scalability is currently limited by transverse mode instabilities (TMI) \cite{zervas2}. The onset of these TMI is in turn influenced by the intensity noise in such fiber lasers \cite{stihler}. These insights have recently put the  noise characteristics of fiber amplifiers into the focus of several investigations \cite{gierschke, li, zhao, buikema}. In previous studies the relative intensity noise was used as a figure of merit, whereas the information on the absolute noise level with respect to the theoretical boundary remains unknown. \newline
\indent The theoretically possible minimum noise in any laser is defined by the quadrature variance of the coherent state \cite{clerk}, the so called shot noise limit (SNL). It is known that phase insensitive amplification processes add noise to the input state and therefore degrade the noise performance of the system \cite{haus, caves}. This is the well known 3$\,$dB limit of optical amplification \cite{bachor}. While quantum noise as well as its impact have been investigated in various phase sensitive and insensitive amplifier systems bevore \cite{kuo, levenson}, it has not yet been experimentally studied in detail for Yb-doped high power fiber amplifiers. An approach to improve the noise characteristics of amplifier systems could be coherent combining of multiple amplifiers. In this case multiple beams are amplified separately and then recombined to generate a high power output \cite{liuCC}. Recently, the noise performance of these systems has been theoretically investigated down to the quantum mechanical noise limit for the first time, showing that it can outperform the output noise variances of a linear amplifier under certain conditions \cite{mueller}. A different attempt for optimizing the noise properties of amplifier systems is to reduce the noise level of its input state - the seed laser - as far as possible. A precise knowledge of the amplifier noise with respect to the SNL might allow to assess the future potential of high power fiber amplifiers for noise sensitive fields of research like high precision optical frequency comb generation and quantum metrology \cite{udem, brandt}.\newline
\indent High power fiber systems usually consist of several amplification stages \cite{beier, guiraud, zhao2}. In this contribution, we examine the first amplification stage of such an Yb-doped fiber amplifier chain. The noise transfer functions of rear amplification stages have been measured, and it has been shown that they can dampen the relative intensity noise \cite{gierschke}. However, the total amount of noise in such systems remained unclear. Thus, this contribution is focused on the amplifier parts, that are expected to have the highest noise, the laser seed and the first amplification stage, and presents  intensity noise measurements of such a system relative to the shot noise limit. This work is organized in the following manner: In section 2 an overview of the experimental setup is given. Section 3 presents and discusses the results, which are brought to a conclusion in section$\,$4.

\section{Experimental setup}

\begin{figure}[h!]
\centering\includegraphics[width=10cm]{figures/Setup.pdf}
\caption{Schematic diagram of the Yb-doped fiber amplifier with quantum noise measurement setup. PM Fiber: Polarization maintaining fiber, EOM: electro-optic modulator, YDF: Ytterbium doped fiber, PC: pump combiner, C: collimator, HWP: half-wave plate, PBS: polarizing beamsplitter, ESA: electronic spectrum analyzer.}
\end{figure}

As shown in Figure 1, the experimental setup consisted of two parts: the examined laser source and the quantum noise measurement setup. The laser seed source was an external cavity diode laser at 1064 nm that was phase-modulated to avoid Stimulated Brillouin Scattering in subsequent parts of the laser \cite{harish}. After phase modulation, the light, with a power of 5 mW, was amplified in a PLMA-10/125 Ytterbium-doped fiber to an output power of 1.89 W, corresponding to a gain of approximately $26\,$dB. The pump light at a wavelength of 976 nm was coupled in counter-propagating direction by a pump combiner. Fiber connectors after the seed laser and the EOM allowed to examine these parts individually. At the end of this all-fiber polarization maintaining setup, a collimator was used to couple the light to the measurement setup. \newline
\indent The latter consisted of a typical self-homodyne measurement setup \cite{bachor} and an electro-optic modulator used for checking the common mode rejection capability of the system. A half-waveplate (HWP) and polarizing beamsplitter (PBS) after the modulator were used to transform the created polarization modulation to amplitude modulation. Another set of HWP and PBS were deployed to be able to attenuate the optical power during the measurements. A final HWP and PBS split the signal. The splitted signals were detected with two identically constructed, home built detectors. The signals of the two detectors were added or subtracted by in-house produced electronics. In the difference signal, the classical noise contribution is canceled out and only quantum noise of the light field and electronic noise of the detectors are visible. In the sum signal, the quantum as well as the classical intensity noise of both signals are measured \cite{bachor}. Measurements were acquired either with an electronic spectrum analyser (ESA) or alternatively using an oscilloscope for the low frequency regime. To align the last set of HWP and PBS, a sinusoidal signal was modulated onto the laser using the EOM and monitored the difference signal with the ESA. By adjusting the waveplate to equalize the power in the two arms, we were able to suppress the amplitude modulation with a typical common mode rejection of $\ge$30$\,$dB.

\section{Experimental results and discussion}

In order to investigate the intensity noise of our Yb-doped fiber amplifier we characterized the three main components (seed laser, phase modulator and amplifier) separately for different frequency regimes ranging from the low kHz to the MHz domain. Furthermore, we investigated two different seed laser sources, to gain more insight on the contribution of the seed laser intensity noise in fiber amplifiers. Moreover, we relate these noise measurements to the carefully calibrated shot noise limit.

\begin{figure}[]
\centering\includegraphics[width=12cm]{figures/Laser1ESA.pdf}
\caption{Laser and amplifier amplitude noise spectrum with respect to shot noise a) in the 0-100kHz regime with RBW and VBW of 10Hz b) in the 0.1-5 MHz regime with RBW and VBW of 1kHz. Figures c) and d) show the linear noise over optical power dependency of the difference signal (shot noise) for 50kHz and 3MHz, respectively. Figures e) and f) contain the respective quadratic excess noise dependencies of the sum signal.}
\end{figure}

\indent Figure 2a) and b) show frequency resolved intensity noise spectra for the frequency regimes from 0-100$\,$kHz (a) and 0.1-5$\,$MHz (b). The strong, frequency dependent noise increase in the range between 0-100$\,$kHz is caused by the frequency dependent gain of our ESA. Therefore all noise contributions should be considered by their relative value above the SNL. The SNL was measured by subtracting the signals of the two individual detectors, while the excess noise is given by the added signals as discussed in section 2. To confirm the quantum origin of the measured SNL, we attenuated both the sum and the difference signal over a range of 15$\,$dB in optical power. The corresponding results are presented in Figures 2c) - f). To generate the attenuation measurements, zero-span measurements containing 400 datapoints at a given frequency (50$\,$kHz and 3$\,$MHz, as shown in Figures 2c)-f)) with a bandwidth of 1$\,$kHz were used to obtain reliable amplitude noise values for a given optical power. Figures 2c) and d) hint a clear linear dependence of the noise on optical power for the difference signal. In contrast, the quadratic noise to optical power dependency in Figures 2e) and f) is created by the sum signal. This provides us with strong evidence for the quantum and classical nature of the respective signals \cite{walls}.  Returning to Figures 2a) and b), the seed laser clearly showed excess noise above the shot noise limit (SNL) of about 9$\,$dB in the kHz domain, decreasing to about 4$\,$dB in the MHz domain.  Next we introduced the SBS preventing phase modulator to our system. It is evident that it did not generate additional intensity noise, as expected when properly aligned. Finally, the amplifier was introduced. Here, dispersion could transfer the phase noise to intensity noise, but this effect is, given the narrow laser linewidth ($140\,$pm), short fiber length ($4\,$m) and low modulation frequency (maximal $525\,$MHz), negligible. The amplifier created approximately 26$\,$dB of amplification in terms of optical power. Therefore the amplifier was attenuated by 26$\,$dB, in order to be compatible with the power requirements of the quantum noise measurement setup. The noise level of the amplifier was significantly higher than of the seed laser, indicating that the amplification process was a major noise contribution. Following \cite{tunnermann}, the relative intensity noise transfer function of the pump is a low-pass, which explains the observed frequency dependency of the noise gain shown in Figure 2a). At 10 kHz, we measured the attenuated amplifier to have approximately 34$\,$dB of intensity noise above the SNL. In the MHz regime, the amplifier excess noise follows the seed noise spectrum introducing about 6$\,$dB of intensity noise after attenuation. For the noise contribution at full amplifier power, the introduced optical attenuation as well as all losses have to be considered. To get a best estimate of the excess noise above SNL for all frequencies, we measured the noise spectra given in Figure 2a) and b) for range of optical powers and fitted the attenuation behavior independently for the four cases of SNL, seed excess noise, modulator excess noise and the attenuated amplifier excess noise to extract the exact values of excess noise above SNL for the latter. The excess noise  for the unattenuated amplifier was calculated by inserting the full optical power into the fit of excess noise above SNL. Exemplary values for the full amplifier excess noise are given in Table 1.

\begin{table}[]
\centering
\begin{tabular}{c c c}
& Seedlaser Excess Noise & Amplifier Excess Noise \\
Frequency &  ($P_\text{opt}$ = 1.3$\,$mW) (above SNL) &  ($P_\text{opt}$ = 1.89$\,$W) (above SNL) \\
\hline
\hline
10$\,$kHz &  11.9$\,$dB & 62.8$\,$dB \\
50$\,$kHz &  7.8$\,$dB & 50.6$\,$dB \\
90$\,$kHz &  8.2$\,$dB & 49.3$\,$dB \\
1$\,$MHz &  4.5$\,$dB & 42.1$\,$dB \\
4$\,$MHz &  4.2$\,$dB & 41.6$\,$dB \\
\hline
\hline
\end{tabular}
\caption{Exemplary values for unattenuated amplifier excess noise given in dB above SNL}
\end{table}
To interpret the values given in Table 1, we need to consider the theoretical noise limit of the amplification process. This would be bounded by using a shot-noise limited seed laser and a perfect, quantum limited linear amplifier. For this case, we expect 29$\,$dB of noise in our system for an optical gain of 26$\,$dB \cite{bachor}. This means, that we exceed the theoretical limit for the linear amplifier by at least a factor of 12$\,$dB in the MHz regime, at 10$\,$kHz we transcend the theoretical limit by more than 33$\,$dB. \\
\indent As previously discussed in section 1, the current power scaling limit in Yb-doped fiber amplifiers is set by transverse mode instabilities \cite{zervas2}. Typically, our high power amplifiers exhibit TMI with frequencies up to 10$\,$kHz \cite{beier}, thus this frequency regime is especially interesting in our analysis. To investigate this domain in more detail, we exchanged the electronic spectrum analyzer (ESA) with an oscilloscope. The recorded time traces of 1$\,$s with a sampling rate of 1$\,$Ms/s were Fourier transformed to obtain the magnitude of the signals' power spectral density (PSD) as shown in Figure 3a). Using two input ports of the oscilloscope allowed us to monitor the sum and difference signals simultaneously. In agreement with the previous ESA measurements, we found a flat seed excess noise spectrum with approximately 10$\,$dB above SNL and no additional contribution from the phase modulator. The large variance of the traces can most likely be attributed to the theoretical resolution bandwidth of 1$\,$Hz and could be potentially smoothed out by introducing an artificial larger bandwidth into the measurement. Similar to the previous measurements, we attenuated the optical power over a range of approximately 10$\,$dB. We summed 100 consecutive points of the trace allowing us a measurement with a bandwidth of 0.1$\,$kHz. Exemplary scaling curves are shown in Figure 3b) and c), again confirming the quantum nature of the difference signal in contrast to the excess noise. The attenuated amplifier noise spectrum is mainly flat at a level of 30 dB above SNL. However, there are two distinct features in the spectrum: a broad peak of about 10 dB between 4 to 5 kHz and a higher peak in the regime below 0.5 kHz. The latter is also visible in the SNL trace and thus most likely attributed to electronic circuit noise of the detectors. The broader peak on the other hand stems from the amplification process but its precise origin remains unclear. Potentially, the pump diodes driving electronics cause this contribution. Following the above procedure to account for the optical attenuation, we get a total amplifier excess noise of 67.5$\,$dB above SNL at 3$\,$kHz, again by far exceeding the theoretical limit of 29$\,$dB. \\

\begin{figure}[]
\centering\includegraphics[width=12cm]{figures/Oszi.pdf}
\caption{Measurement of low frequency intensity noise power spectral density in the regime between 0-10$\,$kHz (a) and respective scaling curves to show the noise over optical power dependency for the difference (shot noise) (b) and sum (c) signal at 3$\,$kHz.}
\end{figure}
\begin{figure}[]
\centering\includegraphics[width=12cm]{figures/LaserVergleich.pdf}
\caption{The diagrams show the comparison of two different seedlasers (laser 1 and laser 2) before (a, c) and after (b, d) propagation through the fiber amplifier in the kHz (a, b) and MHz (c, d) frequency range.}
\end{figure}

\indent As a final step, we exchanged our seed laser to a comparable ECDL with larger amplitude noise, denoted as laser$\,$2, and repeated our previous measurements. The complete spectral overview for the seed laser as well as attenuated amplifier stage is given in Figure 4. Consistent with previous measurements from this section, the phase modulator does not introduce additional amplitude noise and is thus excluded in further analysis. Figure 4a) and c) present the two different seed lasers for the kHz (a) and MHz (c) domain. It is evident that the noise level of laser$\,$2 exceeded the one of laser$\,$1 by more than 10$\,$dB in the kHz regime. The excess noise of laser 2 decreased for larger frequencies ending up with approximately 2$\,$dB above the level of laser$\,$1 at 5$\,$MHz. While the amplifier conserved the observed spectrum of both lasers in the MHz measurement (Figure 4d)), the noise properties in the kHz regime are vastly different. The noise increase towards low frequencies is less pronounced in the amplifier seeded with laser$\,$2, while at the same time the variance of the noise trace was heavily enhanced. A possible explanation for this behavior would be given by the large enhancement of the seednoise contribution, potentially exceeding the low frequency pump noise dominance given in the amplifier seeded with laser$\,$1. This result demonstrates, that the seed laser intensity noise characteristics can have great influence on the excess noise of the fiber amplifier in the kHz as well as in the low MHz frequency regime.
\section{Conclusion}

We characterized the intensity noise of an Yb-doped fiber amplifier with respect to shot noise for a variety of different frequency regimes. We identified a significant noise contribution above the shot-noise limit (SNL) in the kHz as well as in the MHz regime. Special interest was given to the region below 10$\,$kHz, where we found the highest intensity noise contribution ($\ge$60$\,$dB above SNL) in our investigation. Furthermore we showed that in our case, the choice of seed laser and its intensity noise had a significant impact on the final noise performance of the amplifier. It can be concluded, that commonly used amplifiers still have a great potential of noise reduction.

\section*{Funding}
Fraunhofer and Max Planck cooperation program (PowerQuant)

\section*{Acknowledgments}
The authors acknowledge the support of Oliver Bittel and Lothar Meier.

\section{}
\subsection{}
\subsubsection{}


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

Since the seminal work of Bender and Boettcher on the existence of all-real eigenvalue spectra exhibited by non-Hermitian Hamiltonians, work on the so-called Parity-Time (PT)-symmetric structures have been pursued worldwide, with unexpected and interesting results, opening a new and exciting way of interpreting the consequences of crossing to the forbidden complex domain \cite{Bender1998a,Bender1999,Bender2002,Bendermust,Bendermakingsense,Benderbook}.  A PT-symmetric quantum system is characterized by a potential function invariant under the simultaneous symmetry operations of parity inversion (P) and time reversal (T). Considering the condition for a generic Hamiltonian $H$ to be PT-symmetric, that is, $(PT)^{-1}HPT = H$ one finds that this requirement is equivalent as stating that the complex potential function has a real part that is even while the complex part is odd under $\mathbf{r}\rightarrow - \mathbf{r}$. These authors studied the eigenvalues of the PT-symmetric Hamiltonian $H = p^2 - (ix)^\varepsilon$ as a function of the real parameter $\varepsilon$ and found a symmetry breaking point, depending on this parameter. More specifically, when $\varepsilon \geq 2$ the eigenvalues are real, positive and discrete. On the other hand, in the so-called broken region, i.e., $\varepsilon < 2$, the eigenvalues switch from real to complex values, and now appear in complex pairs of the form $a\pm ib$. Thus, even though the Hamiltonian is PT-symmetric for every $\varepsilon$, below the symmetry breaking point, or the exceptional point, the spectrum of the PT-symmetric Hamiltonian becomes complex \cite{Heiss_2004}. More specifically, the complex conjugate eigenvalue pairs represent non-equilibrium systems with exponential amplification and dissipation, therefore they represent open systems \cite{Benderbook,rotter2009}. However, below the PT phase transition, the system behaves like a closed system due to the all-real energy spectrum and by a new definition of the inner product in the Hilbert space \cite{Bender2002}. It was later shown in a series of papers published by Mostafazadeh that PT-symmetric Hamiltonians actually belong to a subclass of a more general set of operators named pseudo-Hermitian Hamiltonians \cite{Mostafazadeh20021,Mostafazadeh20022,Mostafazadeh20023}.

As waves always behave in a similar way, the notion of exceptional points and PT symmetry has revealed itself to be quite useful in various areas of physics, such as electromagnetic invisibility cloaks \cite{zhu2013}, acoustic scattering \cite{zhu2014} and metasurfaces \cite{sounas2015}, to cite a few. Owing to the analogy between the time-dependent Schr\"odinger equation and the scalar paraxial wave equation, derived from Maxwell's relations, an earlier experimental verification of a PT-symmetric physical model was performed within the context of optics, by using waves with a high degree of coherence, i.e., laser beams, interacting with coupled waveguides \cite{Ruter2010,Christodoulides2007,guo2009}. The potential energy $V(x)$ of a one-dimensional quantum system is interpreted as the complex refractive index $n(x)$ whose imaginary part represents regions in space of gain and loss, distributed in a balanced way. Many unexpected optical effects were unraveled by the incidence of monochromatic radiation into PT-symmetric media as for example in ultra-sensitive sensors \cite{Ren2017}. The huge impact that PT-symmetric ideas have printed upon laser science and technology may be evaluated by the number of emerging areas that there has been since. One remarkable example is the idea of a laser-absorber device \cite{longhi2010, chong2011}. A recent comprehensible review on the applicability of PT optics and non-Hermitian photonics is given in \cite{Longhi2018}. Nowadays researchers have also realized that the application of the ideas of PT-symmetry and exceptional points to optical systems has become an important area of research in integrated photonics \cite{ozdemir2019}.

However, there is one particular important aspect of the optical fields that remained untouched by the new concepts of PT-symmetric quantum mechanics. The fluctuations of the incident radiation field and their influence on the spectral density of the field, on scattering, have been completely ignored in all studies, at the best of our knowledge. It is well known that complete spatial or temporal coherence of radiation does not exist and thus, a theory that includes the random nature of the incident field, should be of some importance to any reader that deals with the wave nature of radiation \cite{wolfbook}. In this way, the scattering of partially coherent radiation by a deterministic periodic medium has been addressed in \cite{Dusek1995,Wolf2013} within the formalism of classical coherence theory. Therefore, considering the huge success of PT-symmetry in practical applications in optics, and furthermore, that field fluctuations are ubiquitous in real photonic systems, in this work we describe a first attempt to elucidate some of the intriguing aspects involving the scattering of a partially coherent radiation field, by a PT-symmetric material with localized regions of gain and loss. Earlier treatments of scattering of coherent radiation by PT-symmetric structures are available \cite{Hurwitz2016,miri2016}. To this end, in section II we describe the model of a PT-structure, based on a pair of particles, representing localized gain and loss materials. We then apply classical coherence theory, based on the Born approximation, to obtain the spectral density changes of a partially coherent field propagating in the far zone, due to the scattering by the two point scatterers. Section III is devoted to the detailed analysis of the spectral density obtained in Section II, providing more clarity into the origins of the spectral changes. Finally, in Section IV, we conclude and discuss further developments.

\section{Physical model}

Let us begin by considering the formalism of classical coherence theory and suppose that all statistical quantities are essentially stationary. We wish to determine the far-field spectral density $S^{(\infty)}(\mathbf{r},\omega)$, at position $\mathbf{r}$ and angular frequency $\omega$, of a partially coherent radiation field that has been scattered by a localized non-Hermitian PT-symmetric structure.

Let us then define the cross-spectral power density of the incident wavefield
$\{ u^{(i)}(\mathbf{r},\omega) \}$,
by $W^{(i)}(\mathbf{r}_1,\mathbf{r}_2,\omega) = \langle
u^{(i)}(\mathbf{r}_1,\omega)^{*}u^{(i)}(\mathbf{r}_2,\omega) \rangle$. Here,
the brackets denote the average over a statistical ensemble of monochromatic
realizations of the incident field. Analogously, the cross-spectral density of the scattered wavefield $u^{(s)}(\mathbf{r},\omega)$ is defined as $W^{(s)}(\mathbf{r}_1,\mathbf{r}_2,\omega) = \langle u^{(s)}(\mathbf{r}_1,\omega)^{*}u^{(s)}(\mathbf{r}_2,\omega) \rangle$ and the two are related by the integral \cite{wolfbook}
\begin{multline}\label{wswi}
	W^{(s)}(\mathbf{r}_1,\mathbf{r}_2,\omega) = \int_V \int_V W^{(i)}(\mathbf{r}_1',\mathbf{r}_2',\omega)F^* (\mathbf{r}_1',\omega)
	F(\mathbf{r}_2',
	\omega) \\
	\times G^*(|\mathbf{r}_1-\mathbf{r}_1'|,\omega)
	G(|\mathbf{r}_2-\mathbf{r}_2'|,
	\omega)d^3 r_1' d^3 r_2',
\end{multline}
where $F(\mathbf{r}, \omega)$ is the potential function proportional
to the refractive index of the material and $G$ represents the Green's function of the system. We will consider a localized scatterer, as depicted in Figure \ref{fig1}, described by the function
\begin{equation}\label{Fr}
    F(\mathbf{r},\omega) = \delta(y)\delta(z)[ (\sigma +i\gamma)\delta(x-a) +(\sigma -i\gamma)\delta(x+a)],
\end{equation}
with $\sigma$ and $\gamma$ being real positive parameters for each individual scatterer, which can be $\omega$-dependent. Since $\gamma$ is associated with the imaginary part of the potential, it describes the non-Hermitian properties of the material. Thus, the scatterer located at $x = a$ has intrisinc loss ($+i\gamma$) while the scatterer located at $x = -a$ has intrinsic gain $-i\gamma$ (see Figure \ref{fig1}). One can easily
verify that the above material function satisfies the PT-symmetry condition:
\begin{equation}
(PT)^{-1}F(\mathbf{r},\omega)PT = F(-\mathbf{r},\omega)^* = F(\mathbf{r},\omega).
\end{equation}
The outgoing free-space Green's function is given by
\begin{equation*}
	G(|\mathbf{r}-\mathbf{r}'|,\omega) = \frac{e^{ik|\mathbf{r}-\mathbf{r}'|}}{|\mathbf{r}-\mathbf{r}'|}
\end{equation*}
where $k=\frac{\omega}{c}$ with $c$ being the speed of light in vacuum. In deriving Equation \eqref{wswi}, the first-order Born approximation was used, so that the scattered spectral density now reads \cite{wolfbook}
\begin{equation}
S^{(s)}(\mathbf{r},\omega) = W^{(s)}(\mathbf{r},\mathbf{r},\omega).
\end{equation}
\begin{figure}[hbt]
	\centering
	\includegraphics[width=0.7\linewidth]{fig1.pdf}
	\caption{(Color online) The model consists of two point scatterers separated by distance $2a$, one having loss $+i\gamma$ (in red) at position $x = a$ and the other having gain $-i\gamma$ (in blue) at position $x = -a$. They are situated in the optical path of a partially coherent beam with coherence length $\Delta$. The incident direction of the beam is represented by the unit vector $\hat{s}_0$ and the scattered far-field direction is represented by the unit radial vector $\hat{s}$.}
	\label{fig1}
\end{figure}
Furthermore, as we are interested in the fields in the far-zone region where $k|\mathbf{r}| \rightarrow \infty$, the Green's function becomes
$G \sim \frac{e^{ikr}}{r}e^{-ik\hat{s}\cdot\mathbf{r}'}$, where $\hat{s} = \frac{\mathbf{r}}{|\mathbf{r}|}$ and $r = |\mathbf{r}|$. Assuming these approximations, the far-field spectral density $S^{(\infty)}(\mathbf{r},\omega)$ is finally obtained
\begin{multline}\label{Sinf}
	S^{(\infty)}(\hat{s},\omega) = \frac{S^{(i)}(\omega)}{r^2}\int_V\int_V \mu^{(i)}(\mathbf{r}_1,\mathbf{r}_2,\omega)	F^* (\mathbf{r}_1,\omega) \\
	\times F(\mathbf{r}_2,\omega)
	e^{-ik\hat{s}\cdot(\mathbf{r}_2-\mathbf{r}_1) } d^3 r_1 d^3 r_2,
\end{multline}
where we have made the assumption that the spectral density of the incident wavefield is independent of position and it is represented by $S^{(i)}(\omega) = W^{(i)}(\mathbf{r},\mathbf{r},\omega)$. The function $\mu^{(i)}(\mathbf{r}_1,\mathbf{r}_2,\omega)$ is the normalized cross-spectral density, also known as the spectral degree of coherence. Since our objective is to construct a simple physical model suitable to describe the essential characteristics of the scattering process, we must choose the spectral degree of coherence and the scattering potential in such a way that the integral in \eqref{Sinf} may be analytically solved. With the choice of $F(\mathbf{r},\omega)$ given by Equation \eqref{Fr}, we now assume a partially coherent incident homogeneous
wavefield of infinite width, whose spectral degree of coherence is given by,
\begin{equation}\label{mui}
\mu^{(i)}(\mathbf{r}_1,\mathbf{r}_2,\omega) = e^{-\frac{(\boldsymbol{\rho}_1 - \boldsymbol{\rho}_2)^2}{2\Delta^2}}e^{ik\hat{s}_0\cdot (\mathbf{r}_2 - \mathbf{r}_1)},
\end{equation}
where $\hat{s}_0$ is the incident field direction, $\hat{s}_0 \cdot \boldsymbol{\rho}_{1,2} = 0$ with $\boldsymbol{\rho}$ denoting a transverse position vector and $\Delta$ is the coherence length of the beam which characterizes the distance between a pair of points in the wavefield such that correlation exists. Equation
\eqref{mui} defines the coherence properties of the incident field,
whose transverse points are correlated in a Gaussian fashion: when
$\Delta \rightarrow \infty$, the wavefield is fully spatially
coherent as $|\mu^{(i)}(\mathbf{r}_1,\mathbf{r}_2,\omega)| = 1$, meaning that any point in the beam is correlated with every other. After substituting Equations \eqref{Fr} and \eqref{mui} into Equation \eqref{Sinf} we obtain the far zone
spectral density
\begin{multline}\label{Smain}
	S^{(\infty)}(\hat{s},\omega) = S^{(i)}(\omega)\frac{2(\sigma^2 +
	\gamma^2)}{r^2}
	\left\{ 1 + e^{ -\frac{2a^2}{\Delta^2} }
	\right. \\
	\left.  \times \left[ \left( \frac{\sigma^2-\gamma^2}
	{\sigma^2+\gamma^2}
	\right)
	\cos(2kas_x) + \left(  \frac{2\sigma\gamma}{\sigma^2 +
	\gamma^2}
	\right)\sin(2kas_x) \right] \right\},
\end{multline}
where $s_x = \sin\theta\cos\phi$ is the $x$ component of the radial unit vector $\hat{s}$ and we have assumed that $\hat{s}_0 = \hat{z}$,
i.e., the incident wavefield is propagating in the positive $z$ direction.
Equation \eqref{Smain} is the central equation of this work. It provides
the spectral density of a partially coherent wavefield that has
been scattered by a PT-symmetric material described by the material
function of Equation \eqref{Fr}.

\section{Far-zone non-Hermitian scattered Spectral Density}

To grasp the physics behind Equation \eqref{Smain}, let us consider some special cases. First, when the structure is Hermitian, $\gamma = 0$ so that
\begin{equation}\label{SH}
S^{(\infty)}_{H}(\hat{s},\omega) = S^{(i)}(\omega)\frac{2\sigma^2}{r^2}\left[ 1+e^{-\frac{2a^2}{\Delta^2}}\cos(2kas_x) \right].
\end{equation}
Equation \eqref{SH} describes the spectral density in a passive system with a real refractive index. It is seen from this relation that the scattered spectrum may differ in general from the incident one depending on its particular state of coherence. When $\Delta \gg 2a$, that is, when the two scatterers are localized in a region within the coherence length of the incident field, an interference pattern is obtained. Otherwise, i.e, when $\Delta \le 2a$, that is, when the scatteres are located in a region approximately equal or larger than the coherence length of the incident wavefield, the interference pattern loses its visibility and eventually vanishes, a very well-known result \cite{wolfbook}.

We see from Equation \eqref{SH} that if the distance between the scatterers satisfy $2ka_n = n\pi$ with $n$ odd, there is no radiation scattered into the $\theta = \frac{\pi}{2}$ and $\theta = \frac{3\pi}{2}$ directions if $\Delta \gg 2a$ (high coherence limit). Figure \ref{fig2} illustrates the polar plot for the Hermitian spectral density with $n =$ 1, 3, 5 and 7 in the incident plane defined by $\phi = 0$ ($s_x = \sin\theta$) and $\phi = \pi$ ($s_x = -\sin\theta$). It is clearly seen in these pictures the inhibition of radiation for the particular values of $\theta = \frac{\pi}{2},\frac{3\pi}{2}$. We will be mainly interested in these two directions because it is around these angles that the spectral density changes more drastically. Note that these directions are in the line of the two scatterers, i.e., they are perpendicular to the direction of the incident radiation. So, we conclude that a passive, Hermitian, configuration does not radiate perpendicularly to the incident direction for these particular values of $a_n$. Note also that the number of intensity maxima in the interference pattern is twice the discrete index $n$ and that the spectral density is symmetric under $\theta \rightarrow \pi-\theta$ ($0\leq\theta \leq \pi$).
\begin{figure}[hbt]
	\centering
	\includegraphics[width=1\linewidth]{fig2.eps}
	\caption{Hermitian spectral density $r^2 S_{H}^{(\infty)}(\hat{s},\omega)/2\sigma^2 S^{(i)}(\omega)$ generated by a Hermitian two-point scatterers with separation distance $a_n = \frac{n\pi}{2k}$ with (a) $n = 1$, (b) $n = 3$, (c) $n = 5$ and (d) $n = 7$. When the distance between the sources are given by $a_n$, the system emits no radiation at $\theta = \frac{\pi}{2}$ and $\theta = \frac{3\pi }{2}$ as can be deduced from the general expression \eqref{Smain} in the high coherence limit $\Delta \rightarrow \infty$. }
	\label{fig2}
\end{figure}
Let us consider $2ka = 3\pi$ fixed from now on and investigate the complex extension of the spectral density for $\gamma > 0$. Figure \ref{fig3} shows the non-Hermitian spectral density given by \eqref{Smain} for three values of $\gamma$. The black line in Figure \ref{fig3} is the Hermitian configuration with $\gamma = 0$ which is also shown in part ($b$) of Figure \ref{fig2}. The blue line indicates the non-Hermitian, PT-symmetric, spectral density with $\gamma = 0.5$ and the red line was obtained with $\gamma = 1.0$. All three plots were generated by fixing $\sigma = 1$ and $k\Delta = 10ka$. This implies that the coherence length $L_c = \Delta$ and the distance between scatterers $L_s = 2a$ are in ratio $\frac{L_c}{L_s} = 5$ so that the scatterers are inside the coherence length of the wavefield and therefore the wavefield vibrations in the region around the scatterers are highly correlated. There are two main effects discernible in this non-Hermitian situation. First, as Figure \ref{fig3} clearly indicates, the intensity maxima rotate in such a way that they tend towards the angle $\theta = \frac{\pi}{2}$ as $\gamma$ increases. This is indicated in the figure by black arrows pointing to the direction of rotation. The second effect is the asymmetrical appearance of scattered radiation at angles $\theta = \frac{\pi}{2}$ and $\theta = \frac{3\pi}{2}$ that was inhibited in the Hermitian configuration. We highlight this new radiation direction by a dotted box in Figure \ref{fig3}. The exact angular rotation is easily calculated from the general relation \eqref{Smain}. The spectral density maximum at $\theta = 0$ (when $\gamma = 0$), for example, deviates from $\theta = 0$ according to the relation
\begin{equation}\label{dev}
    \theta(\gamma) = \sin^{-1}\left[ \frac{1}{2ka}\tan^{-1}\left( \frac{2\sigma\gamma}{\sigma^2-\gamma^2} \right) \right]
\end{equation}
which is plotted in Figure \ref{fig4} for three values of $\sigma$. Equation \eqref{dev} is valid for $0\leq \gamma \leq 1$ and one can see from Figure \ref{fig4} that the deviation is almost linear in this interval. Also, larger values of $\sigma$ lead to smaller deviation angles from zero for a fixed $\gamma$. When $\sigma = \gamma$ we have $\theta(\gamma) = \sin^{-1}\left(\frac{\pi}{4ka}\right)$ which depends only on the separation distance between the scatterers.

\begin{figure}[hbt]
	\centering
	\includegraphics[width=0.7\linewidth]{fig3.eps}
	\caption{(Color online) Normalized spectral density $r^2 S^{(\infty)}(\hat{s},\omega)/2S^{(i)}(\omega)$ as a function of $\theta$. The range $\theta > \pi $ simply means that $\phi = \pi$. Black line: Hermitian scattering described by Equation \eqref{SH}. Blue line: Non-Hermitian, PT-symmetric, spectral density calculated using Equation \eqref{Smain} with
	$\gamma = 0.5$. Red line: Non-Hermitian spectral density with $\gamma = 1$, obtained from \eqref{Smain}. In this plot, $ka = \frac{3\pi}{2}$	and $k\Delta = 10ka$, meaning that the scatterer is
	situated in a region where the wavefield vibrations are highly correlated. The black arrows indicate the sense of rotation of the interference fringes as the non-Hermitian parameter $\gamma$ increases. Inside the dotted box, non-negligible radiation appears which are not present in its Hermitian counterpart. The value $\sigma = 1$ was assumed in all these plots.}
	\label{fig3}
\end{figure}

\begin{figure}[hbt]
	\centering
	\includegraphics[width=0.7\linewidth]{fig4.eps}
	\caption{ Angle (in degrees) of the maximum emitted radiation in the direction $\theta = 0$ (when $\gamma = 0$ initially) as a function of the non-Hermitian parameter $\gamma$. The variation of the angular direction is given by Equation \eqref{dev}. In this figure we have used $ka = \frac{3\pi}{2}$, as in Figure \ref{fig3}. The deflection is independent of the coherence length $\Delta$.}
	\label{fig4}
\end{figure}

This asymmetric directional character of the scattered radiation present in non-Hermitian optical configurations has been the subject of intense research since the first papers dealing with beam propagation in PT-symmetric lattices or coupled waveguides \cite{Makris2008,Gao2018,Ruter2010,Zhu2016,shui2018}. Here we see that if the classical coherence of a wavefield is taken into account, one can destroy these important and unusual non-Hermitian properties. Figure \ref{fig5} illustrates this aspect by showing the spectral density in two different configurations. In the first case, the two scatterers are situated in a region where the optical wavefield is highly coherent, $L_c = 5 L_s$, meaning that the field fluctuations around the two scatterers are correlated according to Equation \eqref{mui}. This situation is represented by the continuous line in Figure \ref{fig5}. The other situation assumes that the wavefield fluctuations are not correlated in the region around the scatterers, $L_c = \frac{L_s}{2}$. In this more realistic scenario, the asymmetric character of the emitted radiation is suppressed by the coherence properties of the wavefield and, as a consequence, the material begins to radiate into previously forbidden regions. Therefore, one finds more possibilities to manipulate the character of the emitted radiation by tuning the coherence properties of the incident wavefield, without changing the geometry of the system.

\begin{figure}[hbt]
	\centering
	\includegraphics[width=0.7\linewidth]{fig5.eps}
	\caption{Non-Hermitian spectral density of partially coherent diffracted wavefields. The continuous line represents the configuration where the two scatterers are situated in a region of high coherence, $L_c = 5L_s$. The dotted line is plotted with $L_c = \frac{L_s}{2}$ which indicates that the wavefield around each scatterer is poorly coherent. The coherence properties of the wavefield thus inhibits the asymmetrical directional character typical of non-Hermitian systems.}
	\label{fig5}
\end{figure}

A nontrivial relation becomes evident by comparing the value of the spectral density at two different directions perpendicular to the incident direction. We define the quantity
\begin{equation}\label{beta}
    \beta(\gamma) = \frac{S^{(\infty)}(\theta = \frac{\pi}{2},\phi = 0)}{S^{(\infty)}(\theta = \frac{\pi}{2},\phi = \pi)}
\end{equation}
as a measure of the ratio between the value of the spectral density at $(\theta = \frac{\pi}{2},\phi = 0)$ and $(\theta = \frac{\pi}{2},\phi = \pi)$ as we change the non-Hermitian parameter $\gamma$. The interesting thing is that we obtain two completely different behaviors depending on the distance between the scatterers. If this distance satisfies the condition assumed in the previous analysis, i.e., $2a_n = \frac{n\pi}{k}$, then $\beta(\gamma) = 1$ \textit{independent of} $\gamma$, $\Delta$ and $\sigma$, as can be easily calculated by using the relations \eqref{beta} and \eqref{Smain}. This effect can be seen in Figures \ref{fig3} and \ref{fig5}. A plot of $\beta(\gamma)$ in this situation would give us a horizontal line starting at $\beta = 1$ and this implies that the magnitudes of the emitted spectral intensity perpendicular to the incident direction are equal. However, if $2a_n \neq \frac{n\pi}{k}$ the function $\beta(\gamma)$ exhibits a minimum value at $\gamma = \sigma$. Part (a) of Figure \ref{fig6} shows the plot of $\beta(\gamma)$ for fixed $L_c = 5L_s$ (high coherence situation). If
$2ka_n = n\pi$ the plot is the line $\beta = 1$ and if $2ka_n \neq n\pi$, $\beta$ shows a \textit{minimum} at $\gamma = 1$ (dashed line) and $\gamma = 5$ (dashed-dotted line). The value of $\beta$ at both minima in part (a) are equal: $\beta(\gamma = \sigma = 1) = \beta(\gamma=\sigma=5)$. This
is because the value of $\beta(\gamma)$ at $\gamma = \sigma$ is independent
of $\sigma$:
\begin{equation}\label{betasigmagama}
    \beta(\gamma = \sigma) = \frac{1 + e^{-2\frac{a^2}{\Delta^2}}\sin(2ka)}{1 - e^{-2\frac{a^2}{\Delta^2}}\sin(2ka)}.
\end{equation}
In the limit $\gamma\rightarrow +\infty$,
\begin{equation}
\beta(\gamma) \sim \frac{1+e^{-2\frac{a^2}{\Delta^2}}\left[\cos(2ka)+\frac{2\sigma}{\gamma}\sin(2ka)\right]}{1+e^{-2\frac{a^2}{\Delta^2}}\left[\cos(2ka)-\frac{2\sigma}{\gamma}\sin(2ka)\right]},
\end{equation}
and so it approaches $1$ for fixed $\sigma$ and $a$. Part (b) of Figure \ref{fig6} is considered at $\sigma = 1$ fixed. The objective of this plot is to show the effect of coherence on the function $\beta$. When the coherence length $\Delta$ of the wavefield is larger than the distance between the scatterers, the function $\beta$ has a minimum at $\beta = \sigma$ with the value of $\beta$ given by Equation \eqref{betasigmagama} and so the effect of increasing $L_c$ compared to $L_s$ is to increase the minimum value of $\beta$ until it gets equal to unity at the low coherence regime $L_c \ll L$. A zero value of $\beta(\gamma)$ corresponds to radiation being emitted mainly around $\theta = \frac{3\pi}{2}$, the direction more closely related to the scatterer with gain (located at $x = -a$).

\begin{figure}[hbt]
	\centering
	\includegraphics[width=0.8\linewidth]{fig6.eps}
	\caption{(a) Plot of $\beta(\gamma)$, defined by Equation \eqref{beta}, as a function of $\gamma$ for fixed $L_c = 5L_s$ (high coherence situation). The curves show the situation where $2ka_n = n\pi$ (continuous line), $\sigma = 1$ (dotted line) and $\sigma = 5$ (dashed-dotted line). In the last two cases, the minimum value of $\beta$ is achieved when $\gamma = \sigma$. (b) Plot of $\beta(\gamma)$ for fixed $\sigma = 1$ and $L_c = 5L_s$ (continuous line), $L_c = 0.25L_s$ (dotted line) and $L_c = 0.25L_s$ (dashed-dotted line).}
	\label{fig6}
\end{figure}

To gain a deeper understanding on the interplay between PT-symmetry
and classical correlations during scattering processes, it is useful to
calculate the cross-spectral density function $W^{(s)}(\mathbf{r}_1,\mathbf{r}_2,\omega)$ of the scattered radiation field directly from Equation \eqref{wswi} with the potential function \eqref{Fr} to obtain
\begin{align}\label{wsr1r2}
 W^{(s)}(\mathbf{r}_1,\mathbf{r}_2,\omega) &= \frac{2(\sigma^2+\gamma^2)e^{ik(r_2-r_1)}}{r_1 r_2} \nonumber\\
 &\times \left\{ W_{S} \cos[ka(s_{1x}-s_{2x})] \right.\nonumber\\
 &+ \left. {} \left(\frac{\sigma^2-\gamma^2}{\sigma^2+\gamma^2}\right)W_{C}\cos[ka(s_{1x}+s_{2x})] \right. \nonumber \\
 &+ \left. {} \left(\frac{2\sigma\gamma}{\sigma^2+\gamma^2}\right)W_{C}\sin[ka(s_{1x}+s_{2x})] \right\}
\end{align}
where $s_{jx}$ is the $x$-component of the unit vector $\hat{s}_j = \mathbf{r}_j/|\mathbf{r}_j|$ for $i \in \{ 1,2 \}$ and
\begin{equation}
\begin{aligned}
  W_S &= W^{(i)}(x_1 = \pm a,x_2 = \pm a,\omega), \\
  W_C &= W^{(i)}(x_1 = \pm a,x_2 = \mp a,\omega),
\end{aligned}
\end{equation}
are the cross-spectral densities of the \textit{incident} field
at the positions of the two scatterers with $W_S$ representing self-correlation and $W_C$, cross-correlation. The presence of $\gamma$
in Equation \eqref{wsr1r2} suggests that the two-point correlation
function depend crucially on the gain/loss properties of the material, changing its amplitude and/or phase according to the specific values of $\sigma$, $\gamma$. Thus, one might conclude that a PT-symmetric medium
may be used to manipulate the correlations properties of the
scattered optical field.

Generally speaking, when one is dealing with the Schr\"{o}dinger equation with a PT-symmetric potential satisfying $V(x) = V(-x)^{*}$, the time evolution of a state is crucially dependent on the spectral properties of the Hamiltonian. It is expected that if the Hamiltonian is dependent upon a parameter, $H(q)$, with $q$ a real number, then there exists a critical value $q_c$ such that the eigenvalues of $H$ are real for $q < q_c$ and form complex conjugated pairs for $q > q_c$, as discussed in the introduction section. The time evolution of a physical state $\psi$ is unitary under a $CPT$ inner product and this is only possible in the regime where all eigenvalues are real \cite{Bender2002}. However, in the present scattering model, we have found no indication of any divergences for any value of $\gamma$ and $\sigma$. One might conjecture that these non-unitary effects are not apparent considering
that our results are obtained from a scattering approach. The usual symmetry breaking phenomenon does become mathematically noticeable when the results are obtained within a scattering matrix approach \cite{alu}. Here, we find no abrupt change
in the spectral functions at all when
crossing the symmetry breaking point $\sigma = \gamma$ and we attribute it to the fact that PT-symmetric systems exhibit differences from scattering problems as they are subjected to different constraints. The question on the
existence and behavior of exceptional points in classical optical
coherence theory is still an open question worth
to be pursued. However, we do not attempt to give a satisfactory answer in the present paper. We hope however that this manuscript may provide a first motivation for the initial studies on the interplay of
PT symmetry and classical coherence effects, which are still lacking in the current literature.

\section{Conclusions}

In conclusion, we have found dramatic non-Hermitian changes occurring in the spectral density and in the cross spectral function of a
scattered wavefield by a PT-symmetric localized material, considering the coherence properties of the source. We have obtained explicit analytic expressions for these  spectral properties under the Born approximation and analyzed many limiting cases to grasp the physics of the interplay between PT-symmetry and partial coherence. In summary, the coherence length of the partially coherent incident beam may alter the asymmetric directional character of the emitted radiation, an effect that has not been previously reported. Furthermore, the cross spectral
correlation exhibits a crucial dependence on the gain/loss spatial
distribution. This feature might be very important for application
purposes due to the possibilities it provides to manipulate the
nuisances of fluctuations by non-Hermitian structures.  As in real life, optical fields are produced by partially coherent sources, and considering the ability of non-Hermitian structures to control and manipulate light transport, we believe that the study on the interaction of PT symmetric structures with partially coherent light should reveal itself quite fruitful with many effects waiting to be explored.

\begin{acknowledgments}
We wish to acknowledge the partial financial support of the Brazilian agencies, CNPq and FAPEAL
\end{acknowledgments}