# Dynamic effects of nuclear surface in isoscalar dipole modes

V. I. Abrosimov<EMAIL_ADDRESS>O. I. Davydovska  Institute for Nuclear
Research, Prospect Nauki 47, 03028 Kyiv, Ukraine

###### Abstract

Dynamic surface effects in the isoscalar dipole modes of heavy nuclei are
studied within a semiclassical model based on the solution of the Vlasov
kinetic equation for finite Fermi systems with a moving surface. In order to
clarify the role of dynamic surface effects we have considered an approximate
solution, which takes into account only the part of the variation of the
phase-space distribution function caused by the dynamic surface (the dynamic-
surface approximation). It is shown that the dynamic surface effects have an
essential influence on the features of the isoscalar dipole modes. The
isoscalar dipole strength function has a two-resonance structure already in
the dynamic-surface approximation, and the centroid energies of both the low-
energy resonance and the high-energy resonance are close to corresponding
centroid energies of exact strength function. Calculations of the velocity
field in the dynamic-surface approximation show the vortex character of the
low-energy isoscalar dipole resonance and the compression character of the
high-energy one.

###### keywords:

isoscalar dipole modes , surface effects , kinetic model , strength function ,
velocity field

###### PACS:

21.60.Ev , 24.30.Cz

##  1 Introduction

Experimental studies of the isoscalar dipole nuclear response have shown not
only the high-energy resonance related to the compressibility of nuclear
matter (the compression mode), but also the low-energy resonance (the vortex
mode) [1]-[4]. Various theoretical studies of nuclear isoscalar dipole
excitations have been carried out within quantum approaches of the random
phase approximation (RPA) type, including their relativistic generalization
[5]-[13], as well as within semiclassical approaches based on the phase-space
dynamics [14]-[19]. It is known that isoscalar dipole excitations of finite
many-body systems are related to the center-of -mass motion of the system and
therefore, in a theoretical study of the isoscalar dipole excitations in
nuclei, the motion of the center of mass should be separated from internal
isoscalar dipole excitations. Moreover, the isoscalar dipole excitations are
considered as the second-order effect in the dipole moment because in the
first order they are reduced to the center-of-mass motion of the system.
Taking into account these features of nuclear isoscalar dipole excitations,
one can expect a significant influence of surface effects on the dynamics of
isoscalar dipole excitations, and therefore it is of interest to study this
problem in more detail.

In this paper, we study the effect of dynamic deformation of nuclear surface
on the isoscalar dipole modes by using a translation-invariant model of small
oscillations of finite Fermi systems based on the Vlasov kinetic equation
[20]. This semiclassical model explicitly uses the collective variables
describing nuclear surface oscillations, which makes it possible to separate
the dynamic surface effects in isoscalar dipole modes and is therefore
suitable for studying our problem. We will focus on the variation of the
phase-space distribution function, which is generated by the dipole dynamic
deformation of the nuclear surface. This variation allows us to obtain
information about the dynamic effects of the nuclear surface in isoscalar
dipole excitations, in particular, by using such physical quantities as the
strength function and the velocity field.

In Section II, the description of the dynamic nuclear surface within our
kinetic model is briefly recalled, while in Section III the dynamic-surface
approximation to the study the surface effects in the isoscalar dipole modes
of nuclei is considered. The results of numerical calculations of the strength
function and the velocity field for isoscalar dipole modes in the dynamic-
surface approximation and the comparison with different approximations are
discussed in Section IV. Section V contains the summary. Mathematical details
on the isoscalar dipole velocity field are given in the Appendix.

##  2 Dynamic nuclear surface in kinetic model

In kinetic model of small oscillations of finite Fermi systems with moving
surface [20], a nucleus is treated as a gas of interacting nucleons confined
to a spherical cavity with moving surface. The dipole surface vibrations are
described by a set of collective variables $\delta R_{1M}(t)$ according to the
equation of surface motion

$R(\vartheta,\varphi,t)=R+\sum\limits_{M}{\delta
R_{1M}(t)}Y_{1M}(\vartheta,\varphi).$ (1)

This equation expresses the time-dependent position of the droplet surface
$R(\vartheta,\varphi,t)$, which is a sphere with radius $R$ in equilibrium.
The variation of the phase-space distribution function $\delta n({\bf r},{\bf
p},t)$ related to the isoscalar dipole excitations is determined by the
linearized Vlasov equation

$\frac{\partial}{\partial t}\delta n({\bf r},{\bf p},t)+\frac{{\bf
p}}{m}\frac{\partial}{\partial{\bf r}}[\delta n({\bf r},{\bf
p},t)-\frac{dn_{0}}{d\epsilon}(\delta V({\bf r},t)+V_{ext}({\bf r},t))]=0$ (2)

and satisfies the boundary conditions at the moving surface that ensure a free
moving surface

$\left[{\delta n({\bf r},{\bf p}_{\bot},p_{r},t)\,-\,\delta n({\bf r},{\bf
p}_{\bot},-p_{r},t)}\right]|\,_{r=R}=-2p_{r}\,\frac{dn_{0}}{d\epsilon}\,\frac{\partial}{\partial
t}\,\delta R(\vartheta,\varphi,t)$ (3)

and

$\delta\Pi_{rr}({\bf r},t)|_{r=R}\,=\,0.$ (4)

Here $p_{r}$ is the radial momentum of the particle, and ${\bf
p}_{\bot}=(0,p_{\vartheta},p_{\varphi})$. The variation of the normal
component of the momentum flux tensor $\delta\Pi_{rr}({\bf r},t)$ is defined
as [21]

$\delta\Pi_{rr}({\bf r},t)=\int d{\bf p}\,p_{r}\,v_{r}\,[\delta n({\bf r},{\bf
p},t)-\frac{dn_{0}}{d\epsilon}\,\delta V({\bf r},t)].$ (5)

The equilibrium distribution function $n_{0}({\bf r},{\bf p})$ depends only on
the one-particle energy $\epsilon({\bf r},{\bf p})$ and we use the
approximation of the Thomas-Fermi type:

$\frac{dn_{0}(\epsilon({\bf r},{\bf
p}))}{d\epsilon}\,=\,-\frac{4}{h^{3}}\delta(\epsilon-\epsilon_{F}),$ (6)

where $\epsilon_{F}$ is the Fermi energy. The mirror reflection boundary
condition (3) determines a connection between the variation of the phase-space
distribution function and surface oscillations, while the boundary condition
(4) provides the consistent motion of particles inside the nucleus and the
free surface.

The residual interaction between nucleons is taking into account in the
kinetic equation (2) and in the boundary condition (4) due to the mean-field
variation $\delta V({\bf r},t)$ given by

$\delta V({\bf r},t)=\int d{{\bf r}}^{\prime}d{\bf p}\,v({\bf r},{{\bf
r}}^{\prime})\,\delta n({{\bf r}}^{\prime},{\bf p},t).$ (7)

In order to obtain an explicit solution of the linearized Vlasov equation, we
assume a separable effective interaction of the form

$v({\bf r},{{\bf
r}}^{\prime})=\kappa_{1}\sum\limits_{M}{r\,{r}^{\prime}Y_{1M}(\vartheta,\varphi)\,Y_{1M}^{\ast}({\vartheta}^{\prime},{\varphi}^{\prime})},$
(8)

where $\kappa_{1}$ is the isoscalar dipole interaction strength.

Taking into account that isoscalar dipole excitations of finite Fermi systems
are a second-order effect with respect to the dipole moment, we assume that
collective isoscalar dipole modes are excited by a weak external field of the
type

$V_{ext}({\bf r},t)=\beta\,\delta(t)\,Q^{(3)}(r)\,Y_{10}(\vartheta,\varphi),$
(9)

where $Q^{(3)}(r)=r^{3}$ is the second-order dipole moment, $\delta(t)$ is the
Dirac delta-function in time, and $\beta$ is a parameter that describes the
external field strength.

To find the solution to dynamic equations (1) and (2) with boundary conditions
(3) and (4), it is reasonable to change from phase variables $({\bf r},{\bf
p})$ to new variables $(r,\epsilon,l,\alpha,\beta,\gamma)$. Here, $\epsilon$
is the particle energy, $l$ its angular momentum, and $\alpha,\beta,\gamma$
are the Euler angles, which describe the rotation to the coordinate system
with the $z$-axis directed along the vector $\bf l={\bf r}\times{\bf p}$ and
the $y$-axis directed along the vector ${\bf r}$. The (Fourier transformed in
time) variation of the distribution function can be written in terms of new
variables as the series expansion [22]

$\displaystyle\delta n({\bf r},{\bf p},\omega)$ $\displaystyle=$
$\displaystyle\sum\limits_{MNL}\,[\delta
n_{MN}^{L+}(\epsilon,l,r,\omega)+\delta n_{MN}^{L-}(\epsilon,l,r,\omega)]$
(10) $\displaystyle\times$
$\displaystyle(D_{MN}^{L}(\alpha,\beta,\gamma))^{\ast}\,Y_{LN}(\frac{\pi}{2},\frac{\pi}{2}),$

where $L=1$ and $M=0$ for the variation of the distribution function related
to the isoscalar dipole excitations induced by the external field (9), the
quantities $D_{MN}^{L}(\alpha,\beta,\gamma)$ are the Wigner rotation matrices
and the value of $\omega$ defines the excitation energy of the nucleus
($E=\hbar\omega$). The functions $\delta n_{MN}^{L\pm}$ refer to particles
having positive ($\delta n_{MN}^{L+}$) or negative ($\delta n_{MN}^{L-}$)
components of the radial momentum $p_{r}$.

##  3 Dynamic-surface approximation

Taking into account the boundary condition (3) the solution of the linearized
Vlasov equation (2) for the (Fourier transformed in time) variation of the
phase-space distribution function can be written as

$\delta n({\bf r},{\bf p},\omega)\,=\,\delta n_{stat}({\bf r},{\bf
p},\omega)\,+\,\delta n_{dyn}({\bf r},{\bf p},\omega).$ (11)

Here $\delta n_{stat}({\bf r},{\bf p},\omega)$ is the solution for a system
confined by the static (fixed) surface while the variation $\delta
n_{dyn}({\bf r},{\bf p},\omega)\sim\delta R_{1M}(\omega)$ is due to surface
oscillations. Taking into account the expansion (10), the solution (11) can be
rewritten in terms of new variables in the form

$\delta n_{N}^{\pm}(r,\epsilon,l,\omega)=\,\delta
n_{stat,N}^{\pm}(r,\epsilon,l,\omega)\,+\,\delta
n_{dyn,N}^{\pm}(r,\epsilon,l,\omega).$ (12)

We can find an explicit expression for the distribution function associated
with isoscalar dipole excitations in the presentation (12) using the equation
(106) of paper [20] with L=1 and M=0111Note that, instead of
$\delta{\tilde{f}}_{0N}^{1\pm}(r,\epsilon,l,\omega)$, here the notation
$\delta n_{N}^{\pm}(r,\epsilon,l,\omega)$ is used. and then obtain, in
different approximations, the isoscalar dipole response function, as well as
the local dynamic quantities associated with isoscalar dipole modes, in
particular, the velocity field .

To evaluate the effect of nuclear surface oscillations on the isoscalar dipole
modes we neglect the static-surface contribution $\delta n_{stat}({\bf r},{\bf
p},\omega)$ to the solution $\delta n({\bf r},{\bf p},\omega)$ and consider
the approximation (the dynamic-surface approximation) defined as

$\delta n({\bf r},{\bf p},\omega)\,\approx\,\delta n_{dyn}({\bf r},{\bf
p},\omega),$ (13)

or in terms of new variables

$\delta n_{N}^{\pm}(r,\epsilon,l,\omega)\,\approx\,\delta
n_{dyn,N}^{\pm}(r,\epsilon,l,\omega).$ (14)

We use this approximation in the energy region 10 - 30 MeV where the isoscalar
dipole modes are observed in heavy spherical nuclei [1, 2, 3, 4].

In kinetic model, the isoscalar dipole response function for a system with a
moving surface is defined as [16, 18]

$R(\omega)=\beta^{-1}\int\,d{\bf
r}\,r^{3}\,Y_{1M}^{\ast}(\vartheta,\varphi)\,\delta{\bar{{\varrho}}}({\bf
r},\omega).$ (15)

Here the modified particle density variation $\delta{\bar{{\varrho}}}({\bf
r},\omega)$ is given by

$\delta\mathop{\bar{{\varrho}}}({\bf r},\omega)=\delta\varrho({\bf
r},\omega)+\varrho_{0}\,\delta(r-R)\,\delta R_{1M}(\vartheta,\varphi,\omega),$
(16)

where

$\delta\varrho({\bf r},\omega)=\int\,d{\bf p}\,\delta n({\bf r},{\bf
p},\omega).$ (17)

In our translation-invariant model the strength associated with the center of
mass motion is concentrated at the zero energy. Therefore the response
function (15) can be written in the form

$R(\omega)=R_{intr}(\omega)-R_{c.m.}(\omega).$ (18)

Here the response of the center of mass $R_{c.m}(\omega)$ induced by the
external field (9) is proportional to $1/{{\omega}^{2}}$. It does not excite
the system at ${\omega}\neq 0$. Taking into account the expantion of the
solution given by the expression (11), the internal response function
$R_{intr}(\omega)$ can be written in the form

$R_{intr}(\omega)=R_{stat}(\omega)+R_{dyn}(\omega)-R_{c.m.}(\omega),$ (19)

where $R_{stat}(\omega)$ is the response function of the system with a fixed
surface and the second term $R_{dyn}(\omega)$ is the contribution to the
internal response function generated by the moving surface. Below we will
consider the internal isoscalar dipole strength function that is determined as

$S(E)=-\frac{1}{\pi}ImR_{intr}(E),$ (20)

where $E=\hbar\omega$. By using the expansion (19) we can evaluate the
strength function in different approximations, in particular, in the dynamic-
surface approximation.

Next, we consider the velocity field associated with the isoscalar dipole
modes. This local dynamic quantity describes the spatial distribution of the
average nucleon velocity under collective excitation and gives the information
about the nature of this excitation. In kinetic theory, the velocity field is
determined as

${\bf u}({\bf r},\omega)=(m\rho_{0})^{-1}\int{d{\bf p}\,{\bf p}\,\delta n({\bf
r},{\bf p},\omega)},$ (21)

where $\rho_{0}$ is the nuclear equilibrium density. We will consider the
velocity field in the meridian plane $XZ$ that usually exploited in the
random-phase-approximation (RPA) calculations [23]. In this representation the
particle radius-vector is given by ${\bf r}=(x,y=0,z)$ or ${\bf
r}=(r,\vartheta,\varphi=0)$ in the spherical coordinates and the velocity
field (21) can be written as

${\bf u}(r,\vartheta,\varphi=0,\omega)=u_{r}(r,\vartheta,\omega){\bf
e}_{r}+u_{\vartheta}(r,\vartheta,\omega){\bf e}_{\vartheta},$ (22)

where $u_{r}(r,\vartheta,\omega)$ and $u_{\vartheta}(r,\vartheta,\omega)$ are
the radial and tangential components of the velocity field vector and ${\bf
e}_{r}$, ${\bf e}_{\vartheta}$ are the corresponding unit vectors of the
spherical coordinate system. The cartesian componens
$u_{x}(r,\vartheta,\omega)$ and $u_{z}(r,\vartheta,\omega)$ of the velocity
field vector (22) are given in Appendix.

The expressions for the functions $u_{r}(r,\vartheta,\omega)$ and
$u_{\vartheta}(r,\vartheta,\omega)$ can be found in the form (see Appendix)

$u_{r}(r,\vartheta,\omega)=Y_{10}(\vartheta,0)u_{r}(r,\omega),$ (23)
$u_{\vartheta}(r,\vartheta,\omega)=Y_{11}(\vartheta,0)u_{\vartheta}(r,\omega).$
(24)

Here the dipole spherical harmonics $Y_{1m}(\vartheta,0)$ determine the
angular dependence of the velocity field, while the radial functions
$u_{r}(r,\omega)$ and $u_{\vartheta}(r,\omega)$ can be written as

$u_{r}(r,\omega)=\pi\frac{1}{\rho_{0}}\frac{1}{r^{2}}\int{d\epsilon}\int{dl}l\sum\limits_{N=-1}^{1}{[\delta
n_{N}^{+}(r,\epsilon,l,\omega)-\delta n_{N}^{-}(r,\epsilon,l,\omega)]},$ (25)
$\displaystyle u_{\vartheta}(r,\omega)$ $\displaystyle=$
$\displaystyle-\frac{i}{\sqrt{2}}\pi\frac{1}{\rho_{0}}\frac{1}{r^{2}}\int{d\epsilon}\int{dl}l\frac{p_{\bot}(r,l)}{p_{r}(r,\epsilon,l)\,}$
$\displaystyle\times$ $\displaystyle\sum\limits_{N=-1}^{1}{N[\delta
n_{N}^{+}(r,\epsilon,l,\omega)+\delta n_{N}^{-}(r,\epsilon,l,\omega)]},$

where $p_{r}(r,\epsilon,l)=[2m\epsilon-{(l/r)}^{2}]^{1/2}$ and
$p_{\bot}(r,l)=l/r$ are radial and tangential momentum of the particle,
respectively. By means of the expansion (12), we get the expression of the
velocity field (22), which allows for the study of the velocity field in
various approximations, in particular, in the dynamic-surface approximation.

##  4 Surface effects in strength function and velocity field

The isoscalar dipole strength function (20) for heavy spherical nuclei in the
energy region of the isoscalar dipole resonances (10-30 MeV) is shown in Fig.
1 in different approximations. We consider the response of a sample ‘nucleus’
of A=208 nucleons. Dipole strength functions calculated for other values of A,
corresponding to other medium-heavy spherical nuclei, are qualitatively
similar to the case shown in Fig. 1, in particular, in Ref. [16], see Fig. 4
(solid curves), it was shown that the isoscalar dipole strength functions for
the spherical nuclei ${{}^{208}Pb}$, ${{}^{116}Sn}$, and ${{}^{90}Zn}$ have a
similar form, but shift to higher energies with decreasing of A. The strength
parameter of the isoscalar dipole residual interaction (8), chosen in order to
reproduce the experimental value of the giant monopole resonance energy in
${{}^{208}Pb}$ within our kinetic model [16], is $\kappa_{1}=-7.5\cdot
10^{-3}$ MeV/fm2. The corresponding value of the incompressibility modulus
equals K${}_{A}=$160 MeV. The numerical calculations were carried out using
the following values for the nuclear parameters: r${}_{0}=$ 1.25 fm,
$\epsilon_{F}=$ 30.94 MeV, and m $=$ 1.04 MeV (10-22sec)2/fm2.

Figure 1: Isoscalar dipole strength function in the energy region of the
isoscalar dipole resonances is shown in different approximations: the exact
strength function within our model (solid curve), the strength function in the
dynamic-surface approximation, see Eq. (13), (dot-dashed curve) and evaluated
with static-surface boundary condition (dotted curve).The system contains
A$=$208 nucleons. The experimental centroid energies of low- and high-energy
isoscalar dipole resonances for ${{}^{208}Pb}$ of Ref. [2] are indicated by
arrows.

By comparing the three curves shown in Fig. 1, we can evaluate the role of
dynamic surface effects to the formation of isoscalar dipole resonances. The
solid curve shows the exact internal strength function within our kinetic
model $R_{intr}(\omega)$, which is calculated using the total variation of the
distribution function $\delta n_{stat}({\bf r},{\bf p},\omega)$ induced by the
external field (9), see (19), (11). It can be seen from Fig. 1 that the exact
strength function has a two-resonance structure. The centroid energies of the
low- and high-energy resonances equal 11.5 and 24.3 MeV, respectively, and are
in reasonable agreement with the experimental centroid energies for the
isoscalar dipole resonances in the ${{}^{208}Pb}$ nucleus [2], which are
indicated by arrows in Fig. 1. The results of calculations of the internal
strength function in the dynamic-surface approximation are shown in Fig. 1 by
the dot-dashed curve. We can see that, in the dynamic-surface approximation,
the strength function also has a two-resonance structure and, moreover, the
centroid energies of the low- and high-energy resonances are very close to the
corresponding centroid energies of the exact strength function. Finally, the
dashed curve in Fig.1 shows that the strength distribution in the static-
surface approximation does not have a significant effect on isoscalar dipole
resonances. Thus, both low-energy and high-energy isoscalar dipole resonances
are formed due to dynamic surface effects.

Figure 2: Isovector dipole strength function taken from Ref. [25] is shown in
different approximations: the strength function with taken into account
dynamic surface effects (solid curve), and the one evaluated with static-
surface boundary condition (dashed curve).The system contains A$=$208
nucleons. The experimental centroid energy of isovector giant dipole resonance
for ${{}^{208}Pb}$ of Ref. [24] is indicated by arrow. Figure 3: Velocity
fields in the $XZ$-plane associated with the low-energy resonance of the
isoscalar dipole strength function at the centroid energy 11.5 MeV, see Fig. 1
(solid curve), are shown in different approximations, see Eqs.(12), (13): the
velocity field evaluated using the exact solution, Fig. 3(a), in the dynamic-
surface approximation, Fig. 3(b), and in the static-surface approximation,
Fig. 3(c). The system contains A$=$ 208 nucleons. Figure 4: Velocity fields in
the $XZ$-plane associated with the high-energy resonance of the isoscalar
dipole strength function at the centroid energy of 24.3 MeV, see solid curve
in Fig.1, in the same approximations as in Fig. 3.

It is interesting to compare the effect of dynamic surface on isoscalar dipole
resonances with the effect on the isovector giant dipole resonance. In Fig. 2
taken from Ref. [25], it can be seen that, in contrast to isoscalar dipole
resonances, dynamic surface effects lead only to a shift of the isovector
giant dipole resonance to lower energies , while this resonance is formed in a
Fermi system with a fixed surface.

To evaluate the effect of the dynamic surface on the character of the velocity
fields associated with isoscalar dipole resonances, in Figs. 3 and 4 are shown
in different approximations the velocity fields at centroid energies of the
low- and high-energy resonances of the strength function. In Figs. 3(a) and
4(a), the exact velocity fields within our kinetic model, see Eqs. (11),(21),
are shown. As can be seen, the exact velocity field associated with the low-
energy resonance at the centroid energy 11.5 MeV displays a purely vortex
(toroidal) character, see Fig. 3(a), while the exact velocity field associated
with the high-energy resonance at the centroid energy 24.3 MeV has a
compression character, see Fig. 4(a). These results are in qualitative
agreement with the previous results of relevant RPA calculations [7, 10, 12,
13].

In Figs. 3(b) and 4(b), the results of calculations of the velocity fields for
the isoscalar dipole resonances in the dynamic-surface approximation are
shown. Fig. 3(b) shows that, in the dynamic-surface approximation, the
velocity field associated with the low-energy resonance has essentially a
vortex character similar to the exact velocity field, cf. Figs. 3(b) and 3(a),
but this velocity field also involves a compression inside the nuclear fluid.
On the contrary, in the static-surface approximation, the velocity field for
the low-energy resonance has no vortex motion, see Fig. 3(c), however it has a
visible strengthening inside the nuclear fluid associated with single-particle
excitations, which ensures a purely vortex character of the exact velocity
field. Therefore, the vortex nature of the low-energy isoscalar dipole
resonance is associated with both dynamic surface effects and single-particle
excitations inside the nuclear fluid. In Fig. 4(b) the velocity field
associated with the high-energy isoscalar dipole mode is shown in the dynamic-
surface approximation. It can be seen that, in the dynamic-surface
approximation, this velocity field has a compression character. Moreover, it
looks very similar to the exact velocity field, cf. Figs. 4(b) and 4(a), while
in the static-surface approximation, the velocity field for the high-energy
isoscalar dipole mode is negligible, see Fig. 4(c). Thus, the features of the
velocity field associated with the high-energy isoscalar dipole mode are
mainly determined by dynamic surface effects .

##  5 Summary

We have studied the effect of the nuclear surface oscillations on the
isoscalar dipole resonances in heavy spherical nuclei within the kinetic model
of small vibrations of finite Fermi systems with a moving surface. The
isoscalar dipole strength function and the velocity fields associated with
isoscalar dipole modes were evaluated in different approximations, in
particular, in the dynamic-surface approximation, which takes into account
only the variation of the phase-space distribution function caused by the
moving surface.

It is found that, in the energy region of the isoscalar dipole resonances
(10-30 MeV), the strength distribution is mostly determined by dynamic surface
effects. The isoscalar strength function of the dipole, evaluated in the
dynamical surface approximation, is very similar to the strength function ,
which is exact within our model. Moreover, without taking into account dynamic
effects of the nuclear surface, both the low-energy isoscalar dipole resonance
(the vortex mode) and the high-energy one (the compression mode) cannot be
created.

Dynamic surface effects also determine the character of the velocity fields
associated with the isoscalar dipole resonances. In the dynamic-surface
approximation, the velocity field for low-energy resonance has a vortex
character. At the same time, the velocity field associated with high-energy
resonance has a compression character that is approximately the same as the
velocity field, which is exact within our model.

The isoscalar dipole resonances (vortex and compression modes), which are
observed in heavy spherical nuclei are essentially related to dynamic effects
of the nuclear surface. Moreover, the formation of vortex motion in low-energy
isoscalar dipole resonance is due to both dynamic surface effects and single-
particle excitations inside a nuclear liquid.

For the further investigation of the nature of collective excitations in
nuclei, it would be interesting to study local dynamic quantities associated
with changes in momentum space, in particular, the momentum flux tensor, which
contains information about dissipative processes (viscosity, thermal
conductivity) in a nuclear quantum fluid.

Acknowledgments. This work is supported in part by the budget program ”Support
for the development of priority areas of scientific researches” of National
Academy of Sciences of Ukraine, project No. 0122U000848.

##  Appendix A Isoscalar dipole velocity field

In this Appendix we give some details of the derivation of the analytical
expressions (23) and (24) for the isoscalar dipole velocity field.

Taking the $z$-axis in the direction of the external field, we are interested
in the velocity field (21) in the meridian plane $XZ$. First we find the
cartesian coordinates of the velocity field $u_{x}({\bf r},\omega)$ and
$u_{z}({\bf r},\omega)$:

$u_{x(z)}({\bf r},\omega)=\frac{1}{m\varrho_{0}}\int d{\bf p}p_{x(z)}\delta
n({\bf r},{\bf p},\omega)\,.$ (27)

It is useful to express them through the spherical components $u_{1m}({\bf
r},\omega)$:

$\displaystyle u_{x}({\bf r},\omega)=-\frac{1}{\sqrt{2}}[u_{11}({\bf
r},\omega)-u_{1-1}({\bf r},\omega)],$ $\displaystyle u_{z}({\bf
r},\omega)=u_{10}({\bf r},\omega),$ (28)

that are related to the corresponding components of the momentum in the lab
system,which can be written as

$p_{1m}(\vartheta_{p},\varphi_{p})=pC_{1m}(\vartheta_{p},\varphi_{p})\,,$ (29)

with $p=\sqrt{2m\epsilon}$ and
$C_{1m}(\vartheta_{p},\varphi_{p})=\sqrt{\frac{4\pi}{3}}Y_{1m}(\vartheta_{p},\varphi_{p})$.

Now we make the same rotation of the coordinate system as in the expansion
(10). In the rotated system,
$p^{\prime}_{1m}(\vartheta^{\prime}_{p},\varphi^{\prime}_{p})=pC_{1m}(\vartheta^{\prime}_{p},\varphi^{\prime}_{p})$.
Note that $\vartheta^{\prime}_{p}=\frac{\pi}{2}$, while
$\frac{\pi}{2}\leq\varphi^{\prime}_{p}\leq 3\frac{\pi}{2}$. The spherical
components of the momentum in the rotated system are conveniently expressed in
terms of the cartesian coordinates:

$\displaystyle
p_{11}(\vartheta^{\prime}_{p},\varphi^{\prime}_{p})=-\frac{1}{\sqrt{2}}(p_{x^{\prime}}+ip_{y^{\prime}}),$
$\displaystyle
p_{10}(\vartheta^{\prime}_{p},\varphi^{\prime}_{p})=p_{z^{\prime}},$
$\displaystyle
p_{1-1}(\vartheta^{\prime}_{p},\varphi^{\prime}_{p})=\frac{1}{\sqrt{2}}(p_{x^{\prime}}-ip_{y^{\prime}})\,,$
(30)

which, for given $\epsilon$ and $l$, are given in terms of the magnitude of
the radial $p_{r}$ and tangential $p_{\bot}$ components of the particle
momentum:

$\displaystyle p_{z^{\prime}}=0,$ $\displaystyle p_{y^{\prime}}=\pm
p_{r}=\pm\sqrt{2m\epsilon-\frac{l^{2}}{r^{2}}},$ $\displaystyle
p_{x^{\prime}}=-p_{\bot}=-\frac{l}{r}\,.$ (31)

According to Eq. (2.24) of [26] the spherical components of the particle
momentum in the two reference frames are related by

$p_{1m}(\vartheta_{p},\varphi_{p})=\sum_{n}{\Big{(}}{\cal
D}^{1}_{mn}(\alpha,\beta,\gamma){\Big{)}}^{*}p_{1n}(\vartheta^{\prime}_{p},\varphi^{\prime}_{p})\,.$
(32)

Now we consider the velocity field (21), written as

${\bf u}({\bf r},\omega)=\frac{1}{m\varrho_{0}}\int d{\bf
r}^{\prime\prime}\int d{\bf p}{\bf p}\delta({\bf r}-{\bf
r}^{\prime\prime})\delta n({\bf r}^{\prime\prime},{\bf p},\omega).$ (33)

We change the integration variables from phase variables $({{\bf
r}}^{\prime\prime},{\bf p})$ to new variables
$(r,\epsilon,l,\alpha,\beta,\gamma)$ according to [22]

$\int d{\bf r}^{\prime\prime}\int d{\bf p}=\int dr^{\prime\prime}\int
d\epsilon\int d\frac{l}{v_{r}}\int_{0}^{2\pi}d\int_{0}^{\pi}\sin\beta
d\beta\int_{0}^{2\pi}d\gamma,$ (34)

take the usual multipole expansion of the $\delta$-function

$\delta({\bf r}-{\bf
r}^{\prime\prime})=\frac{\delta(r-r^{\prime\prime})}{r^{2}}\sum_{L^{\prime}M^{\prime}}Y_{L^{\prime}M^{\prime}}(\vartheta,\varphi)Y_{L^{\prime}M^{\prime}}^{*}(\vartheta^{\prime\prime},\varphi^{\prime\prime})\,,$
(35)

and the transformation analogous to Eq. (32)

$Y_{L^{\prime}M^{\prime}}^{*}(\vartheta^{\prime\prime},\varphi^{\prime\prime})=\sum_{N^{\prime}}{\cal
D}^{L^{\prime}}_{M^{\prime}N^{\prime}}(\alpha,\beta,\gamma)Y_{L^{\prime}N^{\prime}}^{*}(\frac{\pi}{2},\frac{\pi}{2})\,.$
(36)

Moreover, we use the expansion (10) for given $L$ and $M=0$.

Then the spherical components of the velocity field (33) are given by:

$\displaystyle u_{1m}^{L}({\bf r},\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{m\varrho_{0}}\sum_{N}\sum_{L^{\prime}M^{\prime}N^{\prime}}Y_{L^{\prime}M^{\prime}}(\vartheta,\varphi)Y_{L^{\prime}N^{\prime}}^{*}(\frac{\pi}{2},\frac{\pi}{2})Y_{LN}(\frac{\pi}{2},\frac{\pi}{2})$
(37) $\displaystyle\times$ $\displaystyle\int
dr^{\prime\prime}\frac{\delta(r-r^{\prime\prime})}{r^{2}}\int d\epsilon\int
dl\frac{l}{v_{r}}\delta n^{L}_{0N}(\epsilon,l,r^{\prime\prime},\omega)$
$\displaystyle\times$
$\displaystyle\sum_{n}p_{1n}(\epsilon,l,r^{\prime\prime})\int_{0}^{2\pi}d\alpha\int_{0}^{\pi}d\beta\sin\beta\int_{0}^{2\pi}d\gamma$
$\displaystyle\times$ $\displaystyle{\cal
D}^{L^{\prime}}_{M^{\prime}N^{\prime}}(\alpha,\beta,\gamma){\Big{(}}{\cal
D}^{1}_{mn}(\alpha,\beta,\gamma)){\Big{)}}^{*}{\Big{(}}{\cal
D}^{L}_{0N}(\alpha,\beta,\gamma){\Big{)}}^{*}\,.$

The angular integrals are given by (see Ref.[26], p. 148)

$\displaystyle\int_{0}^{2\pi}d\alpha\int_{0}^{\pi}d\beta\sin\beta\int_{0}^{2\pi}d\gamma{\cal
D}^{L^{\prime}}_{M^{\prime}N^{\prime}}(\alpha,\beta,\gamma){\Big{(}}{\cal
D}^{1}_{mn}(\alpha,\beta,\gamma){\Big{)}}^{*}{\Big{(}}{\cal
D}^{L}_{0N}(\alpha,\beta,\gamma){\Big{)}}^{*}$
$\displaystyle=(-)^{m-n}\frac{8\pi^{2}}{2L+1}<L^{\prime}1M^{\prime}-m|L0><L^{\prime}1N^{\prime}-n|LN>.$
(38)

For given $L$ and $M=0$, the Clebsh-Gordan coefficients vanishes, unless
$|L-1|\leq L^{\prime}\leq L+1$, $M^{\prime}=m$ (this fixes $M^{\prime}$) and
$N^{\prime}=N+n$. Moreover,
$Y_{L^{\prime}N^{\prime}}^{*}(\frac{\pi}{2},\frac{\pi}{2})$ vanishes unless
$L^{\prime}+N^{\prime}$ is even. The components of the particle momentum
$p_{1n}=0$ for $n=0$, thus $n=\pm 1$. This implies $N^{\prime}=N\pm 1$ and
$L^{\prime}=L\pm 1$ (no term $L^{\prime}=L$). Therefore, for given $L$ and
$M=0$, Eq. (37) gives

$\displaystyle u_{1m}^{L}({\bf r},\omega)$ $\displaystyle=$
$\displaystyle(-)^{m}[<(L-1)1m-m|L0>Y_{L-1m}(\vartheta,\varphi)u_{L,L-1}(r,\omega)$
(39) $\displaystyle+$
$\displaystyle<(L+1)1m-m|L0>Y_{L+1m}(\vartheta,\varphi)u_{L,L+1}(r,\omega)],\,$

with the radial functions

$\displaystyle u_{L,L\pm 1}(r,\omega)$ $\displaystyle=$
$\displaystyle-\frac{1}{m\varrho_{0}}\frac{8\pi^{2}}{2L+1}\sum_{N=-L}^{L}\sum_{n=\pm
1}Y_{LN}(\frac{\pi}{2},\frac{\pi}{2})Y_{L\pm
1N+n}^{*}(\frac{\pi}{2},\frac{\pi}{2})\,$ (40) $\displaystyle\times$
$\displaystyle<(L\pm 1)1(N+n)-n|LN>$ $\displaystyle\times$
$\displaystyle\frac{1}{r^{2}}\int d\epsilon\int dl\frac{l}{v_{r}}\delta
n^{L}_{0N}(\epsilon,l,r,\omega)p_{1n}(\epsilon,l,r)\,.$

For $L=1$ the equations (39) and (40) gives the analytical expressions for the
spherical components of the dipole velocity field. They can be written as

$\displaystyle u_{10}^{1}({\bf
r},\omega)=Y_{00}u_{10}(r,\omega)-\sqrt{\frac{3}{5}}Y_{20}(\vartheta,\varphi)u_{12}(r,\omega),$
$\displaystyle u_{1\pm}^{1}({\bf
r},\omega)={\pm}\sqrt{\frac{3}{10}}Y_{21}(\vartheta,\varphi)u_{12}(r,\omega),$
(41)

with the radial functions

$\displaystyle u_{10}(r,\omega)$ $\displaystyle=$
$\displaystyle-i\sqrt{\frac{1}{3}}\pi\frac{1}{\varrho_{0}}\frac{1}{r^{2}}\int
d\epsilon\int dll$ $\displaystyle\times$
$\displaystyle\sum_{N=-1}^{1}\\{i[\delta
n^{1+}_{0N}(\epsilon,l,r,\omega)-\delta n^{1-}_{0N}(\epsilon,l,r,\omega)]$
$\displaystyle+$ $\displaystyle
N\frac{1}{p_{r}(\epsilon,l,r)}\frac{l}{r}[\delta
n^{1+}_{0N}(\epsilon,l,r,\omega)+\delta n^{1-}_{0N}(\epsilon,l,r,\omega)]\\},$
$\displaystyle u_{12}(r,\omega)$ $\displaystyle=$
$\displaystyle-i\sqrt{\frac{2}{3}}\pi\frac{1}{\varrho_{0}}\frac{1}{r^{2}}\int
d\epsilon\int dll$ $\displaystyle\times$
$\displaystyle\sum_{N=-1}^{1}\\{-i[\delta
n^{1+}_{0N}(\epsilon,l,r,\omega)-\delta n^{1-}_{0N}(\epsilon,l,r,\omega)]$
$\displaystyle+$
$\displaystyle\frac{N}{2}\frac{1}{p_{r}(\epsilon,l,r)}\frac{l}{r}[\delta
n^{1+}_{0N}(\epsilon,l,r,\omega)+\delta n^{1-}_{0N}(\epsilon,l,r,\omega)]\\}.$

Finally, taken into account that in the meridian plane $XZ$ the particle
radius-vector is given by ${\bf r}=(r,\vartheta,\varphi=0)$, and the
expressions for the radial and tangential components of the velocity field in
terms of the the cartesian components are defined as

$u_{r}(r,\vartheta,\omega)=\sin(\vartheta)u_{x}(r,\vartheta,\omega)+\cos(\vartheta)u_{z}(r,\vartheta,\omega),$
(44)
$u_{\vartheta}(r,\vartheta,\omega)=\cos(\vartheta)u_{x}(r,\vartheta,\omega)-\sin(\vartheta)u_{z}(r,\vartheta,\omega),$
(45)

we get the analytical expressions (23) and (24) for the isoscalar dipole
velocity field.

## References

  * [1] H.L. Clark,Y.-W. Lui, and D.H. Youngblood, Phys. Rev. C 63, 031301(R) (2001).
  * [2] M. Uchida, H. Sakaguchi, M. Itoh, M. Yosoi, T. Kawabata, H. Takeda, Y. Yasuda, T. Murakami, T. Ishikawa, T. Taki, N. Tsukahara, S. Terashima, U. Garg, M. Hedden, B. Kharraja, M. Koss, B.K. Nayak, S. Zhu, M. Fujiwara, H. Fujimura, K. Hara, E. Obayashi, H.P. Yoshida,H. Akimune, M.N. Harakeh, and M. Volkerts, Phys. Lett. B 557, 12 (2003).
  * [3] D.H. Youngblood, Y.-W. Lui, H.L. Clark, B.John, Y. Tokimoto, and X.Chen, Phys. Rev. C 69, 034315 (2004).
  * [4] M. Uchida, H. Sakaguchi, M. Itoh, M. Yosoi, T. Kawabata, Y. Yasuda, H. Takeda, T. Murakami, S. Terashima, S. Kishi, U. Garg, P. Boutachkov, M. Hedden, B. Kharraja, M. Koss, B.K. Nayak, S. Zhu, M. Fujiwara, H. Fujimura, H. P. Yoshida, K. Hara, H. Akimune, and M.N. Harakeh, Phys. Rev. C 69, 051301(R) (2004).
  * [5] Nguyen Van Giai and H. Sagawa, Nucl. Phys. A 371, 1 (1981).
  * [6] G. Colo, N. Van Giai, P.F. Bortignon, and M.R. Quaglia, Phys. Lett. B 485, 362 (2000).
  * [7] D. Vretenar, A. Wandelt, and P. Ring, Phys. Lett. B 487, 334 (2000).
  * [8] J.Piekarewicz, Phys. Rev. C 62, 051304(R) (2000).
  * [9] M.L. Gorelik and M.H. Urin, Phys. Rev. C 64, 047301 (2001).
  * [10] D. Vretenar, N. Paar, P. Ring, and T. Niksic, Phys. Rev. C 65, 021301(R) (2002).
  * [11] J. Kvasil, N. Lo Iudice, Ch. Stoyanov, and P. Alexa, J. Phys. G: Nucl. Part. Phys. 29, 753 (2003).
  * [12] A. Repko, P.-G. Reinhard, V.O. Nesterenko, and J. Kvasil, Phys. Rev. C 87, 024305 (2013).
  * [13] P.-G. Reinhard, V.O. Nesterenko, A. Repko, and J. Kvasil, Phys. Rev. C 89, 024321 (2014).
  * [14] E.B. Balbutsev, I.V. Molodtsova, and A.V. Unzhakova, Europhys. Lett. 26, 499 (1994).
  * [15] V.M. Kolomietz and S. Shlomo, Phys. Rev. C 61, 064302 (2000).
  * [16] V.I. Abrosimov, A. Dellafiore, and F. Matera, Nucl. Phys. A 697, 748 -764 (2002).
  * [17] M. Urban, Phys. Rev. C 85, 034322 (2012).
  * [18] V.I. Abrosimov and O.I. Davidovskaja, Ukr. J. Phys. 61, 571 (2016).
  * [19] V.I. Abrosimov and O.I. Davydovska, Ukr. J. Phys. 63, 1043 (2018).
  * [20] V.I. Abrosimov, A. Dellafiore, and F. Matera, Phys. Part. Nucl. 36, 699 (2005).
  * [21] E.M. Lifshitz and L.P. Pitaevsky, Physical Kinetics (Pergamon Press, London, 1979).
  * [22] D. M. Brink, A. Dellafiore, and M. Di Toro, Nucl. Phys. A 456, 205 (1986).
  * [23] I. Hamamoto, H. Sagawa, and X.Z. Zhang, Nucl. Phys. A 648, 203 (1999).
  * [24] A. van der Woude, Electric and Magnetic Giant Resonances in Nuclei (World Scientific, Singapore, 1991), Chap. 2, p. 99.
  * [25] V.I. Abrosimov and O.I. Davidovskaja, Bull. Rus. Acad. Sci. Phys. 67, 1757 (2003).
  * [26] D. M. Brink, and G. R. Satchler, Angular momentum (Clarendon press, Oxford, 1968).


# Optimal Grain Mixing is NP-Complete
Tech Report #MSU-NISL-21-001

Md Asaduzzaman Noor<EMAIL_ADDRESS>
Gianforte School of Computing
Montana State University
Bozeman, MT 59717 Sean Yaw<EMAIL_ADDRESS>
Gianforte School of Computing
Montana State University
Bozeman, MT 59717 Binhai Zhu<EMAIL_ADDRESS>
Gianforte School of Computing
Montana State University
Bozeman, MT 59717 John W. Sheppard<EMAIL_ADDRESS>
Gianforte School of Computing
Montana State University
Bozeman, MT 59717

###### Abstract

Protein content in wheat plays a significant role when determining the price
of wheat production. The Grain mixing problem aims to find the optimal bin
pair combination with an appropriate mixing ratio to load each truck that will
yield a maximum profit when sold to a set of local grain elevators. In this
paper, we presented two complexity proofs for the grain mixing problem and
showed that finding the optimal solutions for the grain mixing problem remains
hard. These proofs follow a reduction from the $3$-dimensional matching
($3$-DM) problem and a more restricted version of the $3$-DM known as planar
$3$-DM problem respectively. The complexity proofs do suggest that the exact
algorithm to find the optimal solution for the grain mixing problem may be
infeasible.

Keywords: Grain Mixing, Precision Agriculture, Computational Complexity

##  1 Introduction

Generally, agriculture and agricultural products are essential in sustaining
lives on the planet. Considerable planning is required in the agricultural
sector To feed the large population on the Earth efficiently. In this paper,
we considered an important component in the wheat distribution profit referred
to as grain mixing (wheat blending). When selling wheat in a local grain
elevator, many variables come into play for determining the price of the
wheat. Among them, protein content plays the most important role which is
affected by several environmental factors such as temperature, soil nitrogen
level, precipitation, etc. Due to these variances, protein content in wheat
changes not only from year to year but also from crop to crop.

After the harvesting season, usually, the farmers store their grain into
multiple grain bins and transport the grain via trucks in batches to sell
wheat in the local grain elevators. Today, the device for tracking protein
content in wheat is available, however, it is expensive. Therefore, small
farmers take their harvest to the nearest elevator and collect the price paid
by them. The grain mixing problem aims to determine the optimal mixing of
different quality wheat (in terms of protein content) to load trucks that
maximize the overall profit when selling wheat to multiple local grain
elevators.

This paper extends our previous work (Noor and Sheppard (2021)) where we
applied and adapted two different evolutionary approaches: genetic algorithms
and differential evolution in order to solve the grain mixing problem. The
experimental results suggested that mixing grain increases the profitability
when selling wheat, and the evolutionary approaches consistently led to higher
overall profit. In this paper, we investigated the complexity proof for the
grain mixing problem. We showed that the grain mixing problem is NP-Complete
following a reduction from the ($3$-DM) and a planar ($3$-DM) problem
respectively. Both of these problems are known to be NP-Complete (Garey and
Johnson (1990); Dyer and Frieze (1986)). The complexity proofs justify the use
of approximation algorithms (such as evolutionary approaches) for getting a
feasible solution.

##  2 Background

There have been several approaches in the literature that attempt to solve the
wheat blending problem (or blending problems in general). A few works in the
literature attempted to solve the decision version of the wheat blending
problem using linear programming (LP). Hayta and Cakmalki (2001) utilized LP
methods capable of predicting the optimal wheat blend ratio for a targeted
final quality to produce a bread-making flour. Haas (2011) used the simplex
algorithm to find the optimum blend that satisfies the customer’s specific
solvent retention capacities (SRC). However, for the grain mixing problem
studied here, the protein cost function is non-linear, and there is no
targeted wheat quality (the protein content of a truck is determined in
runtime). Therefore, it can not be applied directly to the LP models.

Mixed Integer Linear programming (MILP) is often used to solve real-world
blending problems with problem-specific constraints. Although the MILP model
can be used to get exact solutions, it is known to be NP-Hard (Krentel
(1988)). Bilgen and Ozkarahan (2007) proposed an MILP model to optimize the
cost for the wheat supply chain (blending, loading, transportation, and
storage), where the model used a specific blending formula for mixing. Meta-
heuristic approaches are also a popular choice for solving blending
optimization problems. Li et al. (2014) proposed a hybrid evolutionary method
to solve the wheat blending problem in Australia. Their problem closely
relates to ours, however, there are some additional constraints in our
problems based on the US wheat market. The real-world blending problems in the
literature do suggest that an exact solution may not be feasible in many
cases. To the best of our knowledge, there is no complexity proof for the
grain mixing (wheat blending) problem.

##  3 Grain Mixing Problem

###  3.1 A Simple Example

First, we will demonstrate how mixing grain can improve profitability when
selling wheat through a simple example. We collected wheat harvesting data
from a local Montana farmer who tracks the protein content of his wheat. In
Montana, when selling wheat, the price per bushel of grain depends on a range
of protein content. Each elevator has a base protein content range for which a
base price is paid. Then the cost model follows a premium-dockage curve where
if a tuck’s protein content is higher than the base protein range a premium
price is paid based on how far it is from the base level, and the price is
docked if the truck’s protein level is lower than the base protein range. In
many cases, the protein content in a bin is short for reaching a higher price
range. Therefore, mixing it with a high protein content bin when loading a
truck might change the average protein content of the truck to reach a higher
price grade which provides the scope for optimization.

Figure 1 shows an example of the grain mixing problem. There are three bins
with different amounts of bushels and protein content in the example. Table 1
shows the example elevator prices where $[11,12)\%$ shows the base protein
range with $\$4$ price. The price is docked for a lower protein range and
increased for a higher protein range. Farmers who do not track the protein
content of wheat would load trucks separately with grain from each bin and the
price they would get would be $(50*\$3+100*\$4+50*\$6)=\$850$. There is a
fixed capacity for the number of bushels a truck can carry which is $100$ in
the example. However, if they were to track the protein content and mix grains
as shown in the figure; load truck one with $50$ bushels from bin one and $50$
bushels from bin two, and load truck two with $50$ bushels from bin two and
$50$ bushels from bin three, the price they will get would be
$(100*\$4+100*\$6)=\$1000$. Therefore, mixing grain in this scenario increases
the profit by $\$150$.

Figure 1: A simple grain mixing example Table 1: Elevator price for grain mixing example Protein range (%) | Price/Bushel ($)
---|---
$[10,11)$ | $3$
$[11,12)$ | $4$
$[12,13)$ | $6$

Besides protein content, two more variables need to be considered in the cost
model. There is a mixing cost associated when mixing grain to load a truck.
The mixing cost depends on the distance between the two bins that were used to
load a truck. The mixing cost is higher when mixing grain from two bins that
are farther apart. Finally, there is a delivery cost associated with the
distance between the truck site and the elevators.

A key challenge in the grain mixing problem is to find the optimal bin pair
combination and the number of bushels drawn from each bin to load tucks that
will yield maximum profit. In this paper, we showed that finding the optimal
bin pair combination with the mixing ratio is an NP-Complete problem.

###  3.2 Problem Definition

In order to derive the complexity proofs from the $3$-DM problems, we had to
assume that each elevator has a fixed amount of bushels that they will accept.
Although this assumption is not present in the general grain mixing problem,
the proof provides us an idea of the hardness of the grain mixing problem
studied here. Therefore, the grain mixing problem can be defined as follows:

Input:

  * •

$B$: Set of bins with each having a capacity, protein content (%), and
elevator-specific delivery cost. Each pair of bins also has a mixing cost.

  * •

$R$: Set of trucks with each having a capacity.

  * •

$M$: Set of elevators with each having a capacity and protein content cost
function (i.e., $ paid per unit with protein content $p$).

The goal of the Grain Mixing ($GM$) problem is to find an assignment of bins
to trucks to elevators and quantity of grain to transfer such that:

  1. 1.

No bin, truck, or elevator capacities are violated.

  2. 2.

The total profit is maximized.

##  4 Complexity Proofs

The complexity proofs for the grain mixing (GM) problem follow from a
reduction from the standard $3$-DM problem and a more restricted version of
the $3$-DM problem known as planar $3$-DM problem respectively. The standard
$3$-DM problem can be defined as follows:

######  Definition 1 ($3$-DM)

Let $X$, $Y$, $Z$ be finite sets with $|X|=|Y|=|Z|=\alpha$, and let
$T\subseteq X\times Y\times Z$ consist of triples such that $x\in X$, $y\in
Y$, and $z\in Z$. $M\subseteq T$ with $|M|=\alpha$ is a valid 3-Dimensional
Matching (3-DM) if for any two distinct triples $(x_{1},y_{1},z_{1})\in M$,
and $(x_{2},y_{2},z_{2})\in M$, the following holds: $x_{1}\neq x_{2}$,
$y_{1}\neq y_{2}$, and $z_{1}\neq z_{2}$.

Similarly, from the standard $3$-DM problem, the restricted planar $3$-DM can
be defined as follows:

######  Definition 2 (Planar $3$-DM)

Let us consider a standard $3$-DM problem. A bipartite graph can be associated
with the 3-DM instance where one side consists of vertices from $X$, $Y$, and
$Z$ and the other side consists of vertices from triples $T$. There is an edge
between the vertices from the sets to the vertices of the triples if and only
if the set element is in the triple. If the associated graph is planar then it
is referred to as planar 3-DM.

The decision version of the (planar) $3$-DM problem is known to be NP-complete
as shown in one of Karp’s NP-Complete problem lists (Karp (1972)). The
decision problem can be defined as “Given a subset of triples $T$ and an
integer $\alpha$, does there exist a $3$-dimensional matching $M\subseteq T$
where $|M|\geq\alpha$”. Consequently, the optimization problem for the $3$DM
can be defined as “Given a subset of triples $T$, find the $3$-dimensional
matching $M\subseteq T$ that maximizes $|M|$”. As the decision problem is NP-
Complete, that follows that the optimization problem will be NP-Hard.

###### Theorem 3

The grain mixing problem is NP-Complete.

Proof  We provide two proofs for the NP-Completeness of the GM problem. The
first one follows a reduction from the planar $3$-DM problem. First, given an
instance of planar 3-DM problem, we convert it to an instance of GM problem as
follows:

  * •

For each element in set $X$, and set $Y$:

    * –

Create Bin $b_{x_{i}}$ and Bin $b_{y_{z}}$ and assign half unit of bushels,

    * –

Set protein content for $b_{x_{i}}$ to $(p-\epsilon)\%$ and for $b_{y_{z}}$ to
$(p+\epsilon)\%$,

    * –

Set initial mixing cost to $+\infty$ for all of the created bins.

  * •

For each element in set $Z$, create an elevator $m_{z_{k}}$ and define its
cost function to pay $\$0$ initially. Set the capacity of the elevators to
accept a maximum of one unit bushel.

  * •

For each triples $(x_{i},y_{j},z_{k})\in T$, create a truck $t_{i}$ that loads
grain from $b_{x_{i}}$ and $b_{y_{z}}$ and deliver to elevator $m_{z_{k}}$.
The maximum truck capacity is also one unit bushel.

    * –

change the mixing cost between $b_{x_{i}}$ and $b_{y_{z}}$ to $C$,

    * –

Set the delivery cost of truck $t_{i}$ to $C$,

    * –

Update elevator $m_{z_{k}}$ cost function to pay $\$R$ if the protein range of
the truck is in $\left[p,p+2\epsilon\right)$.

Figure 2 shows the reduction from the planar 3-DM to the GM instance. In the
figure, the triples represented in the red box could be considered as an
intermediate warehouse that is connected to the elevator by an edge with cost
$+\infty$ that is used to prevent flowing trucks/bushels from one elevator to
another (as it is not possible for the GM scenario). The triples vertex with
the associated set vertex can be thought of as a truck starting at $b_{x_{i}}$
location and going to $b_{y_{j}}$ location for mixing grain incurring mixing
cost $C$ and from there going to the warehouse with a delivery cost of $C$.
These costs are represented with the edge distance. We then prove the
following claim: The planar 3-DM problem has a matching of size $\alpha$ if
and only if the GM instance has a maximum revenue of $\$R\times\alpha$ and a
minimum cost of $2C\times\alpha$.

Figure 2: Reduction from planar 3-DM to GM instance

If part: If the planar 3-DM has a matching of size $\alpha$ then the GM
instance a maximum revenue of $\$R\alpha$ and a minimum cost of $2C\alpha$.

For each matching in triples $(x_{i},y_{j},z_{k})\in\alpha$, we load trucks
with half unit bushels from bin $b_{x_{i}}$, and half unit from bin
$b_{y_{j}}$ and send to elevator $m_{z_{k}}$. The elevator will pay $\$R$ as
the protein range of the truck would be $\left[p,p+2\epsilon\right)$ and the
total mixing cost and delivery cost would be $2C$. Therefore, for $\alpha$
matching, the total revenue would be $\$R\alpha$ and the total cost would be
$2C\alpha$.

Only If part: If the GM instance a maximum revenue of $\$R\alpha$ and a
minimum cost of $2C\alpha$ then the planar 3-DM has a matching of size
$\alpha$.

First, we argue that no two bins from the same set $X$ or $Y$ will be in the
valid GM solution as the mixing cost between the same set of bins is set to
infinity. Therefore, there will be no profit from mixing from the same set of
bins. Then, the same bin will not appear twice in the valid GM solution as it
will either increase the overall mixing cost or decrease the overall revenue.
For example, each bin has a capacity to hold half unit bushels and if a bin is
delivered to multiple elevators the truck will be partially filled. For a
partially filled truck the revenue would be $<\$R$ and there will be multiple
mixing costs for the same bin which would be $>C$. Therefore, using the same
bin to deliver grain to multiple elevators would either decrease the overall
revenue or increase the overall cost.

Finally, the elevator can accept only one unit of bushels. Therefore, if it
accepts a fully loaded truck (truck capacity is also one unit bushels) then it
will not accept bushels from any other truck. However, two partially loaded
trucks may deliver to the same elevator which will violate it to be a valid
matching. We argue that partially filled trucks would not be in the valid GM
solution as in that case it will either increase the overall delivery cost or
decrease the overall revenue. For example, if a partially loaded truck
delivers to an elevator, the elevator will pay $<\$R$ amount and in that case,
another partially filled truck has to deliver to the same elevator to get an
overall revenue of $\$R$. However, in this case, multiple trucks have to go to
the same elevator which will increase the delivery cost to be $>C$. Therefore,
each truck in the valid GM solution will have a unique bin pair and an
elevator providing a planar 3-DM of size $\alpha$. As GM is obviously in NP,
we have the theorem.

Proof  The second proof follows a reduction from the standard $3$-DM problem.
Given an instance of the $3$-DM problem, we turn it into an instance of the
$GM$ problem as follows:

  1. 1.

Create bin $b_{x_{i}}$ for each $x_{i}\in X$ and create bin $b_{y_{j}}$ for
each $y_{j}\in Y$. Give all bins the capacity of one half unit.

  2. 2.

For each bin $b_{a}$, give it a unique protein content $p_{a}$ between $0$ and
$1$ such that the average protein content of each pair of bins is unique. This
can be done by keeping a list of disallowed values including already assigned
protein values, as well as protein values that when averaged to an already
assigned value will equal some pair’s average value. To select a new protein
content value, a random value that is not on the disallowed value list is
generated.

  3. 3.

Set all mixing costs initially to infinity.

  4. 4.

For each $z_{k}\in Z$, create the elevator $m_{z_{k}}$ and give it a capacity
of one unit. Set all elevator cost functions initially to zero (i.e.,
elevators pay $$0$ per unit, regardless of protein content).

  5. 5.

Create $|T|$ trucks, each having a capacity of one unit.

  6. 6.

Assign the following values:

     * •

$\beta=\min\\{|p_{a}-p_{b}|:a,b\text{ are different bins}\\}$. $\beta>0$ since
the protein content of each bin is unique.

     * •

$\delta=\min\\{|\frac{p_{a}+p_{b}}{2}-\frac{p_{c}+p_{d}}{2}|:\text{ bins
}\\{a,b\\}\neq\\{c,d\\}\\}$. $\delta>0$ since the average protein content for
each pair of bins is unique.

     * •

$\omega=\frac{2\delta}{\beta}>=2$ since
$\delta=|\frac{p_{a}+p_{b}}{2}-\frac{p_{c}+p_{d}}{2}|=|\frac{p_{a}-p_{c}+p_{b}-p_{d}}{2}|\geq\frac{\beta+\beta}{2}=\beta$
implies $\frac{2\delta}{\beta}\geq\frac{2\delta}{\delta}=2$.

  7. 7.

Set delivery cost for each bin to $\omega$.

  8. 8.

For each $t=(x_{i},y_{j},z_{k})\in T$:

    1. (a)

Set mixing cost as $mix(b_{x_{i}},b_{y_{j}})=0$.

    2. (b)

Modify $m_{z_{k}}$’s cost function to pay $(1+\omega)$ per unit for protein
content of $\frac{p_{x_{i}}+p_{y_{j}}}{2}$, in addition to whatever the cost
function already was.

Figure 3 shows the reduction from $3$-DM to $GM$. We then proof the following
claim: The planar 3-DM problem has a matching of size $\alpha$ if and only if
the GM instance has a profit of $\$\alpha$.

Figure 3: Reduction from $3$-DM to the GM instance.

If part:  Suppose the $3$-DM instance has a matching $M$ of size $\alpha$. For
each $t=(x_{i},y_{j},z_{k})$ in $M$, put one half unit from bin $b_{x_{i}}$
and one half unit from bin $b_{y_{j}}$ on truck $t$ and deliver to elevator
$m_{z_{k}}$. Since each $x_{i}$, $y_{j}$, and $z_{k}$ can occur at most once
each in $M$, the capacity constraints of bins $b_{x_{i}}$ and $b_{y_{j}}$,
truck $r_{z_{k}}$, and elevator $m_{z_{k}}$ will not be violated. Also, since
the protein content of truck $t$ is $(p_{x_{i}}+p_{y_{j}})/2$, elevator
$m_{z_{k}}$ will pay $(1+\omega)$ for its one unit. The profit of this
assignment is $\sum_{t\in M}(1+\omega-\omega)=\sum_{t\in M}(1)=\alpha$.

Only If part:  Suppose the $GM$ instance has an assignment whose profit is
$\alpha$. Consider a elevator $m_{z_{k}}$ with profit larger than zero. This
elevator must have received material with a protein content matching some
entry of its cost function. In other words, the protein content of its
received material must equal $(p_{x_{i}}+p_{y_{j}})/2$ for some
$(x_{i},y_{j},z_{k})\in T$. We will show that it is not possible for any bins
other than $b_{x_{i}}$ and $b_{y_{j}}$ to have provided this material to
$m_{z_{k}}$ at a profit. Suppose that this material came from bins
$b^{\prime}_{x_{i}}$ and $b^{\prime}_{y_{j}}$ instead of from $b_{x_{i}}$ and
$b_{y_{j}}$. We first aim to quantify the maximum amount of material that
could be on the truck from $b^{\prime}_{x_{i}}$ and $b^{\prime}_{y_{j}}$ that
achieved the protein content of $(p_{x_{i}}+p_{y_{j}})/2$. Since the average
protein content for each pair of bins is unique, this leads to two cases:

Case 1:
$\frac{p_{x_{i}}+p_{y_{j}}}{2}<\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}$

The maximum amount of material from bins $b^{\prime}_{x_{i}}$ and
$b^{\prime}_{y_{j}}$ that yields a protein content of
$(p_{x_{i}}+p_{y_{j}})/2$ must include the full half unit capacity of the bin
with the lower protein content, say bin $b^{\prime}_{x_{i}}$. Otherwise, the
amount of additional material needed from bin $b^{\prime}_{y_{j}}$ to maintain
the weighted average protein content is less than the amount of material from
left in bin $b^{\prime}_{x_{i}}$. The amount of material needed from bin
$b^{\prime}_{y_{j}}$ to achieve a protein content of $(p_{x_{i}}+p_{y_{j}})/2$
is represented as $q$ in the equation,

$\frac{\frac{1}{2}p^{\prime}_{x_{i}}+qp^{\prime}_{y_{j}}}{\frac{1}{2}+q}=\frac{p_{x_{i}}+p_{y_{j}}}{2}$
(1)

Furthermore, since the difference of average protein content between pairs of
bins is at least $\delta$,

$\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}-\frac{p_{x_{i}}+p_{y_{j}}}{2}\geq\delta$
(2)

Combining these equations gives the bound of,

$\delta\leq\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}-\frac{\frac{1}{2}p^{\prime}_{x_{i}}+qp^{\prime}_{y_{j}}}{\frac{1}{2}+q}$
(3)

Solving for $q$ gives the bound of,

$\displaystyle\frac{\frac{1}{2}p^{\prime}_{x_{i}}+qp^{\prime}_{y_{j}}}{\frac{1}{2}+q}$
$\displaystyle\leq\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}-\delta$ (4)
$\displaystyle\frac{1}{2}p^{\prime}_{x_{i}}+qp^{\prime}_{y_{j}}$
$\displaystyle\leq\frac{1}{2}\left(\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}-\delta\right)+q\left(\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}-\delta\right)$
(5) $\displaystyle
qp^{\prime}_{y_{j}}-q\left(\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}-\delta\right)$
$\displaystyle\leq\frac{1}{2}\left(\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}-\delta-p^{\prime}_{x_{i}}\right)$
(6) $\displaystyle
q\left(\frac{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}}{2}+\delta\right)$
$\displaystyle\leq\frac{1}{2}\left(\frac{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}}{2}-\delta\right)$
(7) $\displaystyle q$
$\displaystyle\leq\frac{\frac{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}}{2}-\delta}{2\left(\frac{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}}{2}+\delta\right)}$
(8)

The total amount of material on the truck delivered to elevator $m_{z_{k}}$
can be calculated as,

$\displaystyle q+\frac{1}{2}$
$\displaystyle\leq\frac{\frac{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}}{2}-\delta}{2\left(\frac{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}}{2}+\delta\right)}+\frac{1}{2}$
(9)
$\displaystyle=\frac{\frac{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}}{2}-\delta}{2\left(\frac{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}}{2}+\delta\right)}+\frac{\frac{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}}{2}+\delta}{2\left(\frac{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}}{2}+\delta\right)}$
(10)
$\displaystyle=\frac{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}}{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}+2\delta}$
(11)
$\displaystyle=\left(\frac{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}+2\delta}{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}}\right)^{-1}$
(12)
$\displaystyle=\left(1+\frac{2\delta}{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}}\right)^{-1}$
(13) $\displaystyle\leq\left(1+\frac{2\delta}{\beta}\right)^{-1}$ (14)
$\displaystyle=\left(1+\omega\right)^{-1}=\frac{1}{1+\omega}$ (15)

So, the truck delivers at most $1/(1+\omega)$ units to elevator $m_{z_{k}}$.
Since the elevator will pay $(1+\omega)$ per unit and the delivery cost is
$\omega$, the profit incurred is at most,

$\frac{1}{1+\omega}(1+\omega)-\omega=1-\omega<=-1,\text{ since $\omega\geq
2$.}$ (16)

Thus, $m_{z_{k}}$ cannot profitably receive material with protein content
$(p_{x_{i}}+p_{y_{j}})/2$ from any bins other than $b_{x_{i}}$ and
$b_{y_{j}}$, which corresponds to an entry $(x_{i},y_{j},z_{k})$ in $M$.

Case 2:
$\frac{p_{x_{i}}+p_{y_{j}}}{2}>\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}$

The maximum amount of material from bins $b^{\prime}_{x_{i}}$ and
$b^{\prime}_{y_{j}}$ that yields a protein content of
$(p_{x_{i}}+p_{y_{j}})/2$ must include the full half unit capacity of the bin
with the higher protein content, say bin $b^{\prime}_{y_{j}}$. The amount of
material needed from bin $b^{\prime}_{x_{i}}$ to achieve a protein content of
$(p_{x_{i}}+p_{y_{j}})/2$ is represented as $q$ in the equation,

$\frac{qp^{\prime}_{x_{i}}+\frac{1}{2}p^{\prime}_{y_{j}}}{\frac{1}{2}+q}=\frac{p_{x_{i}}+p_{y_{j}}}{2}$
(17)

Furthermore, since the difference of average protein content between pairs of
bins is at least $\delta$,

$\frac{p_{x_{i}}+p_{y_{j}}}{2}-\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}\geq\delta$
(18)

Combining these equations gives the bound of,

$\delta\leq\frac{qp^{\prime}_{x_{i}}+\frac{1}{2}p^{\prime}_{y_{j}}}{\frac{1}{2}+q}-\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}$
(19)

Solving for $q$ gives the bound of,

$\displaystyle\delta+\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}$
$\displaystyle\leq\frac{qp^{\prime}_{x_{i}}+\frac{1}{2}p^{\prime}_{y_{j}}}{\frac{1}{2}+q}$
(20)
$\displaystyle\frac{1}{2}\left(\delta+\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}\right)+q\left(\delta+\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}\right)$
$\displaystyle\leq qp^{\prime}_{x_{i}}+\frac{1}{2}p^{\prime}_{y_{j}}$ (21)
$\displaystyle
q\left(\delta+\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}-p^{\prime}_{x_{i}}\right)$
$\displaystyle\leq\frac{1}{2}\left(p^{\prime}_{y_{j}}-\delta-\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}\right)$
(22) $\displaystyle q$
$\displaystyle\leq\frac{p^{\prime}_{y_{j}}-\delta-\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}}{2\left(\delta+\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}-p^{\prime}_{x_{i}}\right)}$
(23)

The total amount of material on the truck delivered to elevator $m_{z_{k}}$
can be calculated as,

$\displaystyle q+\frac{1}{2}$
$\displaystyle\leq\frac{p^{\prime}_{y_{j}}-\delta-\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}}{2\left(\delta+\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}-p^{\prime}_{x_{i}}\right)}+\frac{1}{2}$
(24)
$\displaystyle=\frac{p^{\prime}_{y_{j}}-\delta-\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}}{2\left(\delta+\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}-p^{\prime}_{x_{i}}\right)}+\frac{\delta+\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}-p^{\prime}_{x_{i}}}{2\left(\delta+\frac{p^{\prime}_{x_{i}}+p^{\prime}_{y_{j}}}{2}-p^{\prime}_{x_{i}}\right)}$
(25)
$\displaystyle=\frac{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}}{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}+2\delta}$
(26)
$\displaystyle=\left(\frac{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}+2\delta}{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}}\right)^{-1}$
(27)
$\displaystyle=\left(1+\frac{2\delta}{p^{\prime}_{y_{j}}-p^{\prime}_{x_{i}}}\right)^{-1}$
(28) $\displaystyle\leq\left(1+\frac{2\delta}{\beta}\right)^{-1}$ (29)
$\displaystyle=\left(1+\omega\right)^{-1}=\frac{1}{1+\omega}$ (30)

So, the truck delivers at most $1/(1+\omega)$ units to elevator $m_{z_{k}}$.
Since the elevator will pay $(1+\omega)$ per unit and the delivery cost is
$\omega$, the profit incurred is at most,

$\frac{1}{1+\omega}(1+\omega)-\omega=1-\omega<=-1,\text{ since $\omega\geq
2$.}$ (31)

Thus, $m_{z_{k}}$ cannot profitably receive material with protein content
$(p_{x_{i}}+p_{y_{j}})/2$ from any bins other than $b_{x_{i}}$ and
$b_{y_{j}}$, which corresponds to an entry $(x_{i},y_{j},z_{k})$ in $M$.

In addition to pairs of bins not being able to profitably supply material at a
protein content other than their average content, more than two bins can never
be put on the same truck, since the mixing costs of bins from the same $X$ or
$Y$ set is infinite. Therefore, every profitable elevator $m_{z_{k}}$ must
have received material only from the two bins that formed an entry
$(x_{i},y_{j},z_{k})$ in $M$. Let $M$ be the set of all $(x_{i},y_{j},z_{k})$
from each elevator with a profit larger than zero. We now must argue that
there are no repeated $x$, $y$, or $z$ values in the set $M$:

  * •

The profitability of elevator $z$ is calculated as
$((1+\omega)q_{z}-\omega\gamma)$, where $q_{z}$ is the amount of material
elevator $z$ receives and $\gamma\in\mathbb{Z}_{\geq 0}$ is the number of
trucks that delivered to it. Since each elevator can accept up to one unit
total, $q_{z}\leq 1$. Thus, the profitability of an elevator is at most
$(1+\omega-\omega\gamma)$, which is less than zero for any value of $\gamma$
larger than one since $\omega$ is at least two. So, a single elevator cannot
profitably receive from multiple trucks (i.e. pairs of bins) and there will be
no repeated $z$ values in $M$.

  * •

A elevator that has been delivered a full unit must have come from the full
capacity of a single pair of bins. This means that bins that deliver to an
elevator receiving a full unit cannot be used to deliver to multiple elevators
since they all have a capacity of $\frac{1}{2}$. Thus, if all elevators are
delivered a full unit, there are no repeated $x$ or $y$ values in the set $M$.
However, it is possible for elevators to be profitable without receiving a
full unit, so submaximal elevators (i.e., profitable elevators that receive
less than one unit) can exist in valid $GM$ solutions.

Submaximal elevators are a problem since they can lead to elevators sharing
bins, which would make their corresponding $x$ and $y$ values appear multiple
times in $M$ and prevent it from being a matching. We proceed by first showing
that submaximal elevators must share a bin with at most one other submaximal
elevator. We then show that we can turn a $GM$ solution with submaximal
elevators into a valid matching by selecting either of the submaximal
elevators in each pair, and that doing so will only increase the profit of the
solution and value of the matching.

If an elevator $z$ is submaximal, both of its bins did not deliver their full
$\frac{1}{2}$ capacity to that elevator, since the protein content would be
invalid if a single bin delivered its full capacity and the other did not. If
both bins of a submaximal elevator $z$ deliver less than their $\frac{1}{2}$
capacity, then the bins cannot have been used to deliver to two other
profitable elevators $z^{\prime}$ and $z^{\prime\prime}$. To be profitable,
elevator $z$ must receive more than $\frac{2}{3}$ of a unit, since $\omega\geq
2$ and $((1+\omega)q_{z}-\omega)\leq 0$ for all $q_{z}\leq\frac{2}{3}$. So,
the most capacity that can remain in each bin of a profitable submaximal
elevator that drew evenly from each bin is $\frac{1}{6}-\epsilon$. Thus, if
$t^{\prime}$ and $t^{\prime\prime}$ draw the maximum capacity from their non
shared bins, their delivered capacities are
$\frac{1}{2}+\frac{1}{6}-\epsilon<\frac{2}{3}$. Therefore, if both bins of a
submaximal elevator $t$ deliver less than their $\frac{1}{2}$ capacity, then
at most one of those bins can be used by another profitable submaximal
elevator. If that bin was not delivered to another elevator, then $z$ can
easily be made maximal by delivering all of its bins’ capacity to it. If that
bin did deliver to another profitable elevator $z^{\prime}$, that elevator
must also be submaximal. The profit of the solution can be increased by making
either of the submaximal elevators maximal and the other zero (for a combined
profit of one), since the combined profit of the two submaximal elevators
sharing a bin is at most $\frac{1}{2}$:

$\displaystyle(1+\omega)q_{z}-\omega+(1+\omega)q_{z^{\prime}}-\omega$
$\displaystyle=(1+\omega)(q_{z}+q_{z^{\prime}})-2\omega$
$\displaystyle\leq(1+\omega)\left(\frac{3}{2}\right)-2\omega\text{, since
$q_{z}+q_{z^{\prime}}\leq\frac{3}{2}$}$
$\displaystyle=\frac{3}{2}-\frac{1}{2}\omega\leq\frac{1}{2}\text{, since
$\omega\geq 2$}$

This means that since submaximal elevators share bins with at most one other
submaximal elevator, submaximal elevators are uniquely paired together. Since
turning paired submaximal elevators into a single maximal elevator increases
the profit, we can ensure that there are no repeated $x$ or $y$ values in $M$
with a profit $\geq\alpha$. Therefore, $M$ is a matching of the $3$-DM
instance with size $\geq\alpha$.

##  5 Conclusion

In this paper, we investigated the computational complexity of the grain
mixing problem. We presented two proofs following a reduction from the $3$-DM
and a planar $3$-DM problem respectively to show that finding the optimal
solutions for the grain mixing problem is NP-Complete. This paper extends our
previous work (Noor and Sheppard (2021)) where we used evolutionary approaches
to find a quality solution for the grain mixing problem in a feasible time.
The proofs suggest that brute-force methods to find the optimal solution may
be infeasible and justify the use of approximation algorithms (such as
evolutionary approaches) to find a quality solution.

## References

  * Bilgen and Ozkarahan (2007) Bilge Bilgen and Irem Ozkarahan.  A mixed-integer linear programming model for bulk grain blending and shipping.  _International Journal of Production Economics_ , 107(2):555 – 571, 2007.
  * Dyer and Frieze (1986) M E Dyer and A M Frieze.  Planar 3dm is NP-complete.  _J. Algorithms_ , 7(2):174–184, 1986.
  * Garey and Johnson (1990) Michael R. Garey and David S. Johnson.  _Computers and Intractability; A Guide to the Theory of NP-Completeness_.  W. H. Freeman and Co., 1990.
  * Haas (2011) Nikolas Haas.  _Optimizing Wheat Blends for Customer Value Creation: A Special Case of Solvent Retention Capacity_.  MS Thesis, Kansas State University, USA, 2011.
  * Hayta and Cakmalki (2001) Mehmet Hayta and Unsal Cakmalki.  Optimization of wheat blending to produce breadmaking flour.  _Journal of Food Process Engineering_ , 24:179 – 192, 08 2001.
  * Karp (1972) Richard M. Karp.  _Reducibility among Combinatorial Problems_ , pages 85–103.  Springer US, Boston, MA, 1972.
  * Krentel (1988) Mark W. Krentel.  The complexity of optimization problems.  _Journal of Computer and System Sciences_ , 36(3):490–509, 1988.
  * Li et al. (2014) Xiang Li, Mohammad reza Bonyadi, Zbigniew Michalewicz, and Luigi Barone.  A hybrid evolutionary algorithm for wheat blending problem.  _TheScientificWorldJournal_ , 2014:967254, 02 2014.
  * Noor and Sheppard (2021) Md Asaduzzaman Noor and John W. Sheppard.  Evolutionary grain-mixing to improve profitability in farming winter wheat.  In _Applications of Evolutionary Computation_ , pages 113–129. Springer International Publishing, 2021.