\section{Introduction: some concepts of factorization}

Factorization is the fulcrum
of the collider system, which, realizing the famous quote of Archimedes, is designed to lift new worlds from the vacuum, with
high-energy colliding protons providing the lever arm.  It is central to the energy frontier program, entering into nearly every attempt to achieve precision within
the Standard Model, and in recognizing signals of new physics \cite{Heinrich:2020ybq,Caola:2022ayt}.   As our searches become more sophisticated,
they become more sensitive to the dynamics of the final state.   Such sensitivity, however, may violate the
assumptions necessary to prove factorization.   The purpose of this little essay is to present one person's viewpoint on
the basic concepts that have gone into proofs of  collinear factorization at hadron colliders,  the
fundamental results, and some of their strengths and limitations,
to suggest where we can be confident, and where more work remains to be done.  The presentation is
a mixture of intuitive and technical arguments, hopefully in a conversational style.  It is thus by no means intended
to provide a final word, or to be rigorous, but rather to lead to further conversation as the field moves forward into the next era of
collider physics.

 Factorization is the quantum field theoretic realization of the parton model \cite{Feynman:1969ej,Bjorken:1969ja}, in which large momentum transfer cross sections in hadron-hadron (or lepton-hadron) collisions are the result of  the short-distance collisions of constituents within those  hadrons, each Lorentz-contracted and time dilated in the frame of the other.  The word ``factorization" in this context refers to our ability to treat the momentum distributions of the partons as incoherent from the probabilities of their short-distance interactions, if and only to the extent that they collide at relative velocities approaching the speed of light.  In this limit, before the collision each hadron is in a superposition of quantum states that is completely independent of the other hadron.   At high enough energies, all hadrons are ``pristine" at the moment of collision, in the same mixture of states since baryogenesis.  The basic concept here is {\it causality}, in the sense that signals from one hadron cannot reach the other before the moment of collision.   This implies universality of the kinematic distributions of constituents in fractional momenta, the parton distributions.  In collinear factorization, the hard scattering probability function is computed with the momenta of the colliding partons taken lightlike and in the directions of the colliding hadrons.

 The collision takes place over a short time scale, creating new particles and/or patterns of momentum flow that evolve into the observed final state.   Interactions that take place  long after the short-distance process involve smaller momentum transfers, which for suitably defined observables do not alter the overall flow of momentum, or excite quantum numbers in a manner that can neutralize those created the hard scattering.   When summed with sufficient inclusivity, perturbative singularities associated with this evolution cancel, because they do not affect the total probability that defines the set of states that are summed.   This is an expression of the {\it unitarity} of the underlying theory.   Here,  the identifiable class of final states is typically characterized by decay products of transient, heavy degrees of freedom and/or high-$p_T$ leptons or jets.     These may be Standard Model or not, and final states may be visible or dark, so long as we can recognize their presence, directly or indirectly.

In summary, the central  concepts of factorization are causality, unitarity and universality.    Most commonly, universality refers to parton distributions, but in a closely related sense universality also applies to hard scattering functions. For parton distributions, it applies between different processes, say the proton distributions in deep-inelastic scattering for an electron-proton initial state, and in Drell-Yan, or jet production production for hadron-hadron scattering.    For the short-distance hard-scattering factors, we assume a universality among perturbative calculations with different long-distance completions of the theory.  For example, the production of a photon plus gluon from a quark pair is the same in real QCD as in the four-dimensional limit of dimensionally-continued QCD.    The universality of parton distributions allows us to transfer information from deep-inelastic scattering to hadron-hadron scattering and vice-versa.   The universality of short-distance factors allows us to compute hard scattering for incoming partons rather than model hadrons, clearly a much more efficient procedure.

In principle, as we go to higher orders in hard scatterings, simultaneously fitting more precisely determined parton distributions within processes we confidently identify as Standard Model at short distances, we can achieve more precision in the comparison of high-statistics data to our calculations, with an eye out for discrepancies.   At the same time, factorization allows us to calculate with the same precision the ways that new forces and fields would manifest themselves in these same data.  In the best of possible worlds, we may see behavior that is qualitatively different from the Standard Model predictions (asymmetries, etc.)   In all this, uncertainties in data, and in parton distributions are subjects of increasing interest.

Collinear factorization is well-illustrated by
the inclusive cross section for the production of a Drell-Yan pair or more generally, a heavy electroweak boson of
mass $Q$,
\begin{eqnarray}
Q^4 \frac{d\sigma^{pp\to \gamma^*(Q)X}}{dQ^2 }\,=\, \sum_{ab}\int dx_a \,f_{a/p}(x_a) \int dx_b \,f_{b/p}(x_b)
\,\omega^{ab\to \gamma^*(Q)X}(x_aP_A,x_bP_B,Q)\,,
\label{eq:dy-fac}
\end{eqnarray}
with the functions $f_{i/p}(x_i)$ the parton distributions, and $\omega^{ab\to \gamma^*(Q)X}(x_aP_A,x_bP_B,Q)$ the perturbatively
calculable short-distance functions for partons $a$ and $b$ to produce an off-shell photon of mass $Q$.
Similarly, jet cross sections take the generic collinear-factorized form
\begin{eqnarray}
S\, \sigma^{pp\to \{j_i\}X}\,=\, \sum_{ab}\int dx_a \,f_a(x_a) \int dx_b \,f_b(x_b)
\,\omega^{ab\to \{c_i(p_i)\}X} (x_aP_A,x_bP_B,\{p_i,\delta_i\})\, .
\label{eq:jet-fac}
\end{eqnarray}
where $\{p_i,\delta_i\}$ labels the momenta of the observed jets and the paramters $\delta_i$ that define them, jet cone radii, for example.  (Notice that these parameters generally imply limitations on the additional particles labelled by $X$.)
In both cases, the short-distance, perturbative function, $\omega$, is normalized to be dimensionless.   We have suppressed
the factorization scale, which for this discussion we assume is of the order of the scale of renormalization.
Clear generalizations exist to multiple electroweak boson production followed by decay, and to an unlimited set of
``new physics" scenarios, as long as these involve heavy particles and short-distance interactions.

\section{What's involved in factorization proofs}

As suggested above, the basic ideas that go into arguments for factorization are causality and unitarity, with universality as a consequence.   In informal language, the formation of initial states takes place too early to affect the hard scattering, and subsequent interactions take place to late to undo it.   Also as hinted above, the applicability of these ideas very much depends on what we measure in the final states, and how we combine them.

 First a word about what is meant by {\it proof} in this context.   All proofs I know of assume, as noted above, that hard-scattering functions are universal, and that they can be identified in the scattering of elementary quanta, generally in a theory that is infrared-regulated dimensionally or by other means.   In that theory, long-distance dynamics is to be identified through its (regulated) infrared singularities, generically poles in $4-D$, with $D$ the number of dimensions.  When a cross section is factorized in this theory, we can identify the short distance function by a suitable definition and separation of parton distributions, which contain all infrared singularities.   In this procedure, our ability of factorize depends on our ability to identify singular long-distance, infrared (IR) behavior in the massless limit.   Once identified, we investigate whether in each relevant region of phase and loop momentum space, the integrand is either consistent with factorization, or is part of a nonfactorizing class of contributions that cancels without double counting.  In all these arguments, a key tool is that physical infrared singularities of all types are logarithmic in gauge theory cross sections.  As a result, a single power suppression in a sum over logarithmically-divergent integrals leads to an integral that is finite in perturbation theory.   Effective theory treatments organize sets of regions where logarithmic singularities arise from the outset.  The consistency of any such choice is verified in the process of matching to the full theory at each order of perturbation theory at relevant scales.   The success of this matching at all orders requires in general arguments that go beyond the effective theory itself, and are part of its justification.

 For renormalizable, non-gauge theories near four dimensions, like $\phi^4$ and Yukawa, or $\phi^3$ near six dimensions, IR  singularities are all ``collinear", corresponding to splittings and recombinations of descendants of the two incoming partons, and the task just described is by now well understood.   The treatment of collinear singularities in gauge theories, it is worth noting, also requires the treatment of unphysical, longitudinally-polarized vector particles that share the momenta of active partons, at least in covariant gauges.   Many of the early all-orders discussions of factorization avoided this complication by working in axial or other physical gauges.   A dedicated discussion of the role (and factorization) of unphysical collinear gluons is given, for example, in Ref.\ \cite{Collins:1989gx}.    In effective theory treatments, unphysical gluon polarizations can be absorbed into the definitions of the physical fields \cite{Stewart:2009yx,Becher:2014oda}.

 For gauge theories, an additional distinguishing feature is soft singularities.  In essense, soft gluons can be exchanged between constituents of the incoming hadrons, and such exchanges appear to have the potential to violate causality, and hence factorization.\footnote{Actually, the violation of factorization would be the least of our problems.}    Roughly speaking, it's as if the quarks in proton A could attract the antiquarks in proton B before they collide, even if they approach each other at the speed of light.  They do not do so, because until the collision the fields involved are gauge-equivalent to zero field \cite{Basu:1984ba}, but this feature is obscured in  perturbation theory, which incorporates sums over physical and unphysical degrees of freedom.    Factorization proofs rely on our ability to organize the phases associated with the external charged particles.   The challenge is to disentangle the unphysical component of these phases from the physical gluon-mediated scattering, and to discover when and how the phases cancel in cross sections.   At this stage, we emphasize again that all these arguments depend crucially on treating the momenta of incoming particles, the partons and hadrons, as lightlike.

Early all-orders discussions of factorization in hadron-hadron scattering \cite{Sachrajda:1978ja,Amati:1978wx,Amati:1978by,Ellis:1978ty,Libby:1978qf,Efremov:1978xm} dealt most exhaustively with collinear singularities, which control evolution, also providing a variety of arguments for the cancellation of soft singularities through virtual-real emission cancellation.
This work was inspired in part by a pioneering low-order demonstration of  soft-gluon cancellation in Ref.\ \cite{Politzer:1977fi} and the older $\phi^4$ fragmentation analysis in Ref.\ \cite{Mueller:1974yp}.  In what follows, I will give an informal account of the soft gluon interactions in hadronic hard scattering, which set QCD and other gauge theories apart, and were the primary concern of a second round of factorization studies
\cite{Collins:1983ju,Bodwin:1984hc,Collins:1985ue,Collins:1988ig}.\footnote{For a relatively brief review of the basic arguments, see \cite{Collins:1998ps}, and for a much more complete discussion, Ref.\ \cite{Collins:2011zzd}.}

 \subsection{Soft gluons and phases}

Associated with each parton that arrives at the hard scattering from the past is a nonabelian phase, sensitive only to that component of the gauge field that contracts with the momentum of the parton.
Such interactions can be organized into a phase by the application of an appropriate gauge transformation. Likewise, for the particles that emerge from the hard scattering, their experience of the color fields of other emerging particles is entirely pure gauge, again involving only one component of the gauge field.   These phases differ in an important way.    This difference illustrates the fundamental relationship between causality and analyticity in quantum field theory, and is a direct consequence of the kinematics of the process.
The interactions of fast-moving partons over a semi-infinite  time scale, $\lambda$, take one of the two general forms, in terms of a ``current", labelled here as simply $J$,
\begin{eqnarray}
\tilde J^\mu_+(k\cdot \beta) \ &=&\ \beta^\mu \int_0^\infty d\lambda\ e^{ i \lambda k\cdot \beta} J(\lambda\beta)\, ,
\nonumber\\[2mm]
\tilde J^\mu_-(k\cdot \beta) \ &=&\ \beta^\mu \int_{-\infty}^0 d\lambda\ e^{ i \lambda k\cdot \beta} J(\lambda\beta)\, ,
\end{eqnarray}
where $\beta^\mu=p^\mu/p^0$ defines the light-like velocity of the parton.
The vector $\tilde J^\mu_+(k\cdot\beta)$ is analytic in the upper half-plane in $k\cdot \beta$, simply because, for Im($k\cdot\beta$) and $\lambda$ both positive,
the integral that defines $\tilde J^\mu_+$ becomes exponentially suppressed for large $\lambda$.  Similarly,
$\tilde J^\mu_-(k\cdot\beta)$ is analytic in the lower half-plane of $k\cdot\beta$.

An essential feature of factorization proofs is the wholesale continuation of momentum space contour integrations away
from singularities, into the regions of analyticity just identified.  Along the deformed contours, the ``eikonal approximation", $p\cdot k\gg k^2$, holds for
generic soft momentum $k^\mu$ and high-energy parton momentum $p^\mu$.   This condition, in turn, implies that
in a subdiagram consisting of collinear lines whose momenta are all of order $p$, all components of momentum
$k^\mu$ can be neglected in the analysis of singularities, {\it except} for $\beta\cdot k$.  This is because a single factor like $k^2/(\beta\cdot k)$ renders
a logarithmically divergent integral finite. For lines that terminate at the hard scattering
(``active" incoming partons) or originate at the hard scattering (produced high-$p_T$ partons), we can argue that no additional singularities prevent this deformation \cite{Collins:1989gx} and that all potentially ``acausal" effects
are organized into unitary Wilson line operators (nonabelian phase operators), coming in from the past to the origin, with lightlike
velocities $\beta_i$, or going out to the future with lightlike velocities $\beta_j$,
\begin{eqnarray}
\Phi_{\beta_i,{\rm in}}^{(R_i)}{}(x)
\ &=& \
{\rm P}\ \exp \left[ -ig\int_{-\infty}^0\ d\lambda\,
\beta_i\cdot
A^{(R_i)}(\lambda\beta_i+x)\, \right]\, ,
\nonumber\\[2mm]
\Phi_{\beta_j,{\rm out}}^{(R_j)}{}(x)
\ &=& \
{\rm P}\ \exp \left[ -ig\int^\infty_0\ d\lambda\,
\beta_j\cdot
A^{(R_j)}(\lambda\beta_j+x)\, \right]\, .
\end{eqnarray}
Here $R$ is the color representation of the parton in question,  and P denotes path ordering.    For elastic S-matrix elements in wide-angle scattering, these operators
fully describe the interactions of soft gluons, and generate all associated phases \cite{Sen:1982bt,Feige:2014wja,Ma:2019hjq}.  For hard-scattering cross sections, however, the situation is more complicated.

\begin{figure}[t]
\centering
\vspace{10mm}
\centerline{\figscale{Amplevel_2}{8cm}}
\caption{Representation of components of a hard-scattering amplitude in terms of two incoming and a representative outgoing jet in a typical hard-scattering process.  The labelled gluons attach to (A) an incoming active parton, (B) a parton emerging from the hard scattering and (C) a spectator parton.  Double lines represent the couplings of the unphysically-polarized collinear gluons referred to in the text.  \label{fig:amplevel}}
\end{figure}

 Referring to Fig.\ \ref{fig:amplevel},
 integrals over soft gluon attachments to spectator lines can extend from negative to positive infinity on the same line, and hence inolve both $\Phi_{\rm in}$ and $\Phi_{\rm out}$ operators.
 This leads to a mismatch in analyticity properties between gluons in the initial state and those in the final state, and this is a physical effect.
 This distinction cannot be removed by a gauge transformation, and corresponds physically to gluons that arrive at the time of the hard scattering.
 Such gluons are the origin of the ``Glauber", or ``Coulomb", gluon problem in hadron-hadron scattering \cite{Collins:1981ta,Bodwin:1981fv}.   In essence, Glauber gluons occur in regions of momentum
 space where integration contours are pinched between singularities associated with incoming and outgoing phase operators, hence for which the eikonal approximation fails.
If such a gluon also attaches to the oppositely moving incoming line, we have an essential obstacle to factorization.   Corrections of this type only begin at three loops at the finite level,
with logarithmic enhancemennts starting at four loops.

It is at this stage that the construction of the cross section comes in.   For appropriately defined cross sections, final-state phase operators (the $\Phi_{\rm out}$s)
cancel in the sum over final states.

The purest example of an outgoing phase is in  deep-inelastic scattering, where it is generated by soft gluon interactions with the outgoing quark (and its collinear radiation).   This is the origin of the matrix element definition of parton distributions, which can be written for quark field $q_a(x)$ with flavor index $a$, in the conventional form \cite{Collins:1981uw}
\begin{eqnarray}
f_{a/h}(x,\mu)\ =\ \frac{1}{4\pi}\, \int_{-\infty}^\infty d\lambda e^{-i\lambda n\cdot (xp)}\ \left \langle h(p) \left |
\bar q_a (\lambda n)\, \Phi^\dagger{}_{n,\rm out}(\lambda n)\, \frac{n\cdot \gamma}{2}  \Phi_{n, {\rm out}}(0)   q_a(0) \right |h(p) \right \rangle\, .
\label{eq:pdf}
\nonumber\\
 \end{eqnarray}
 Here, soft singularities cancel in the final state because of the unitarity of the Wilson line operators, which are both oriented in direction $n^\mu$, conventionally chosen
 as a lightlike direction opposite to the lightlike momentum $p$ of incoming hadron $h$.   Argument $\mu$ is the factorization scale, the scale associated with renormalization of the composite operator that defines the matrix element \cite{Collins:1981uw}, which becomes ultraviolet divergent
 because its fields are all separated by lightlike distances ($n^2=0$).   In effect, the scale $\mu$ cuts off the transverse momenta of radiation.   The transverse momentum can
 actually be fixed in this matrix element.   This possibility is realized in transverse momentum distributions and factorization, with many applications in deep-inelastic scattering and Drell-Yan proceses \cite{Rogers:2015sqa}.

 We can extend collinear factorization in deeply inelastic scattering far beyond the total cross section, to sums designed to identify jet and related final states, because collinear
 and soft singularities cancel in the corresponding partial sums over states.   This is a consequence of the  ``local" unitarity \cite{Sterman:1993hfp,Sterman:1995fz,Capatti:2022tit}  of quantum field theories.
 Unitarity at fixed values of spatial momenta is the basis of the infrared safety of jet cross sections in $\rm e^+e^-$ annihilation \cite{Sterman:1978bi,Sterman:1978bj,Sterman:1979uw},  and is closely related to the ideas of effective field theory, and to efforts to reformulate calculations
 entirely in four dimensions \cite{TorresBobadilla:2020ekr}.   For example, an interaction Hamiltonian that excites and absorbs partons in only two directions is not rotationally-invariant, yet it is hermitian, and hence generates a unitary theory.   As such, summing over its final states gives, by the optical theorem, a forward scattering amplitude that is free of
 infrared  singularities, either soft or collinear.

 \subsection{Unitarity, final states and soft gluons}

Moving on to hadron-hadron scattering, the essential observation is that in inclusive cross sections, final state interactions, including phases cancel in the sum over final states.   This will be the case not only
for particles emerging from the hard scattering at large $p_T$, but also for spectator particles that are collinear to the incoming directions.   It is this cancellation that eliminates pinches in the Glauber regions.

The ``local" unitarity of QCD, mentioned above, enables the cancellation for sectors in momentum space that describe jet cross sections.
The long-distance phases  that remain are generated only by initial state operators, $\Phi_{\rm in}$.   Given what we've seen, this is somewhat surprising at first, because it is just the opposite to deep-inelastic scattering, built, as in Eq.\ (\ref{eq:pdf}),  on $\Phi_{\rm out}$ operators. Nevertheless,  universal factorization into parton distribution functions, as naturally defined in deep-inelastic scattering, is still possible so long as soft singularities cancel in the hadron-hadron cross section in question.
 This cancellation was referred to as the transition from ``weak" to ``strong" factorization in the context of Drell-Yan cross sections,  \cite{Collins:1983ju,Bodwin:1984hc,Collins:1985ue,Collins:1988ig}.  Reference \cite{Aybat:2008ct}  discussed explicitly how soft radiation sensitive to the difference between the incoming and outgoing Wilson lines cancels.

 The technical consequence of this cancellation is that the interactions of all ``acausal" initial-state gluons connected to the
set of partons collinear to incoming lines A and B become coherent, and are summarized by two Wilson lines of the type $\Phi_{\rm in}$.
In a fully inclusive cross section, we expect all soft singularities associated with these operators to cancel, enabling a universal factorization.   This is the case of Drell-Yan processes.   It is natural to ask, however,
about less inclusive cross sections.   How inclusive is inclusive enough, and what sums over final states will cancel outgoing nonabelian phases up to a momentum scale at which perturbation theory can be applied?

Here is one way of thinking about these questions.  We can characterize the length scales associated with the hard scattering by the energy scales of heavy or high-$p_T$ particles in the final state.  This alone, however, does not ensure that final state phases cancel.   Rather, a  ``figure of merit" corresponding to this cancellation is  the  coordinate space distance between the hard scattering in the amplitude and in the complex conjugate.   At zero separation, outgoing Wilson lines cancel identically.

The relative ``displacement" of the hard scattering is determined by the range of momenta integrated in the final state, by the usual complementarity of position and momentum measurements.   For hard scale $Q$, if we integrate over a range of momenta in the final state characterized by scales $\delta Q$, with $\delta$ a fraction such that $\alpha_s(Q) \ln (1/\delta) < 1$, we anticipate that the hard scatterings are localized on the scale $1/(\delta Q)$.   Such scales can be used to define jet and related cross sections.   In these cases, the final-state phases $\Phi_{\rm out}$ corresponding to outgoing partons will be matched with their hermitian conjugates at scales that are moderately close to $1/(\delta Q)$, and we can expect their soft gluons to cancel up to the scale $\delta Q$, which is consistent with the definition of the cross section.   As $\delta$ decreases, of course, the cross section becomes more and more sensitive to soft radiation.   At whatever scale the non cancellation sets in, it is no longer possible to absorb radiation at that energy scale (or higher) into cancelling nonabelian initial state phases, and radiation is no longer incoherent between the descendants of the incoming particles.   At scales beyond this one, factorization is ``broken", in the sense found for amplitudes with forward collinear lines in the final state \cite{Catani:2011st}.    Perhaps more accurately, what we have is collinear factorization at an anomalously low factorization scale, determined by the parameters of the set of final states, rather than the hard scattering itself.   Beyond this scale, in fact, virtual or real radiation coupling to fragments of the colliding hadrons, the spectators, can influence measured quantities.  As noted above, such corrections, which are consistent with collinear factorization for inclusive cross sections, begin at three loops.    Important examples of cross sections that illustrate the limitations of collinear (and transverse) momentum factorization in cases like these have been identified, for example, in Refs.\ \cite{Forshaw:2012bi,Rogers:2010dm,Zeng:2015iba,Rothstein:2016bsq,Schwartz:2018obd}.

Given these general considerations, I'll go on in the next section to give a viewpoint on some of the commonly-discussed cross sections, and the arguments I'm aware of that determine their factorization status.

\section{Where we stand}

At the risk of introducing evanescent terminology, I'd like to distinguish between ``deep" factorization, at the scale of the hard scattering, and ``shallow" factorization, in which the factorization scale is set by relatively small (if still perturbative) energy scales that are built into the final state sum in question.    It is important to remember, though, that when such scales decrease, corresponding to a less inclusive cross section, with fewer final states, the cross section must decrease.   To see this decrease, it will be necessary to generalize the analysis of Sudakov resummation, which can be carried out for fixed-angle exclusive amplitudes  \cite{Sen:1982bt,Feige:2014wja,Ma:2019hjq}.   In time, we will learn to love cross sections with shallow factorization, which are in some sense intermediate between inclusive cross sections and exclusive amplitudes.   They encode the coherence of color flow in the perturbative regime, in manners that we have not yet fully understood.

\subsection{Processes with deep factorization}

The arguments above provide, I hope, sufficient framework to interpret at least the claims in the literature of what factorization arguments purport to show.   In these terms, a list of processes for which factorization scales are similar to the corresponding scale of the hard scattering includes electroweak annihilation (Drell-Yan, electroweak boson and Higgs production), and jet production in a wide range of kinematics.

\begin{itemize}

\item  {\it Electroweak annihilation (fully inclusive in hadronic final states). }  These cross sections, were discussed explicitly in the ``second round" of factorization papers, from the mid 1980's \cite{Bodwin:1984hc,Collins:1985ue,Collins:1988ig}.   In this case, the cancellation of final states is complete, a claim that is made explicit by the use of time-ordered or light cone-ordered perturbation theory.   Later work \cite{Aybat:2008ct} explored more in detail the manner in which soft radiation associated with the initial state phases cancels, leading to classical ``universality" for the parton distributions.   Similar results can hold in multiple hard scatterings, where, for example, two electroweak bosons are created at separate partonic interactions in a single collision \cite{Diehl:2015bca}.

\item  {\it Inclusive jet production}.  Here, we imagine sets of states with jets defined through parameters $\epsilon$ with $\alpha_s(Q) \ln (1/\epsilon) \sim 1$.   For a review, see \cite{Sapeta:2015gee}.   The discussion of jet cross sections began early \cite{Sachrajda:1978ja}, and in Ref.\ \cite{Libby:1978qf} the local unitarity mentioned above was used to show that collinear factorization for jet cross sections of this sort reduces to the Drell-Yan case:  soft quanta that connect active or spectator partons to lines in the outgoing jets cancel in the sum over jet final states (assuming that soft radiation is summed inclusively up to a finite fraction of the jet energy scales).   This leaves only soft radiation connecting the two incoming jets, just as in Drell-Yan.

\item  {\it Single-particle inclusive cross sections}, with inclusivity defined similarly to the jet cross sections defined above, provides another example.   Arguments for the use of universal fragmentation functions  were given in Ref.\ \cite{Nayak:2005rt}, very much relying on the concepts described above.

\item   {\it Multi-particle inclusive cross sections for massive quarks}:   As an example, consider final states with a pair of top quarks,.  We may observe the top momentum with arbitrary accuracy, so long as we average over the momenta of the anti-top \cite{Mitov:2012gt}.
For cross sections with the hard scale far above the mass of a heavy quark, additional considerations enter  \cite{Beneke:2016jpx,Caola:2020xup}.
Interestingly, the same cancellations will occur if we measure the transverse momentum of both elements of the top pair, if the pair recoils against sufficiently hard gluonic radiation, over which we average.   This is possible, of course, only at an order high enough that such radiation is present in the final state.     Thus, we expect the cancellation to fail at low orders, and it does \cite{Mitov:2012gt,Collins:2007nk,Collins:2007jp}.   Nevertheless, in the production of very heavy pairs, the cross section is actually dominated by final states with substantial radiation.   In this case, the nonradiative portion of the cross section is Sudakov suppressed.   In this limit, then, it might be possible to compute the pair transverse momentum self-consistently.

\end{itemize}

\subsection{Shallow factorization}

The generic examples of shallow factorization are nonglobal sums over final states \cite{Dasgupta:2001sh}.   In these cross sections, we impose an energy cutoff in one angular region at a scale much lower than that of unobserved ``jets and subjets" in another region.   In general, at high orders these lead to ``super-leading" logarithmic enhancements \cite{Forshaw:2012bi,Becher:2021zkk}, sensitive to both soft and collinear limits, and the factorization scale decreases to the measure of energy flow, rather than the underlying hard scale.   Again, in this case, we can think of these logs as due to the non cancellation of outgoing nonabelian phases for energy scales larger than the cutoff imposed to define the cross section.

Going to the limit, if we ``empty out" a section of phase space (say, a complete rapidity gap), we break the cancellation of final-state interactions, including the matching of nonabelian phases in the amplitude and complex conjugate in hadron-hadron scattering.  Nevertheless, in deep inelastic scattering a variant factorization is possible \cite{Collins:1997sr}, which goes beyond the present discussion.

\section{Conclusions}

In this short review, I've tried to summarize some of the elements of the existing literature on collinear factorization, concentrating on basic intuition and the role of soft radiation.
I have left to the side many topics, including  transverse momentum factorization \cite{Rogers:2015sqa}, factorization for generalized parton distributions \cite{Belitsky:2005qn},  the uses of factorization in resummation, higher twist factorization, and its application in nuclear targets and collisions.
Even purely collinear factorization is a vast subject, and in closing I can only emphasize that even within this narrow focus on a familiar but consequential topic, much important work has surely been neglected.
Renewed interest in this subject, partly inspired by the high energy behavior of the spontaeneously broken sectors of the Standard Model \cite{Bauer:2017isx,Han:2020uid,Chien:2018ohd,Baumgart:2018ntv}, will hopefully lead
to further insights in how to use examples of ``shallow" factorization, to study the transition between short- and long-distance dynamics in gauge theories.

\subsection*{Acknoledgements}  I am grateful to Werner Vogelsang for helpful comments, and to Fabio Maltoni, Shufang Shu and Jesse Thaler for their encouragment of the idea of preparing a short discussion of the status of collinear factorization in connection with the Snowmass process.   This work was supported in part by the National Science Foundation, award PHY 1915093.


\section{Introduction}

	The aim of this paper is to show the existence of at least three weak solutions for the following class of nonlocal quasilinear elliptic systems in $\mathbb{R}^N$

	\begin{eqnarray}
		\label{s1.1}
		\begin{cases}
			-M_i\Big(\mathcal{B}_i(u_i)\Big)
			\textrm{div}\,\Big(  \mathcal{A}_{1i}(\nabla u_i)+ w_i(x)\mathcal{A}_{2i}(u_i)\Big)=\lambda F_{u_i}(x,u) & \text{in }\mathbb{R}^N, \\
			\qquad u_i \in W_{w_i}^{1,p_i(x)}(\mathbb{R}^N)\cap W_{w_i}^{1,\gamma_i(x)}(\mathbb{R}^N) ;
		\end{cases}
	\end{eqnarray}
	for $1\leq i \leq n$  ($n\in \mathbb{N} $), where $N\geq 2$,  $\lambda$ is a positive real parameter and the functions $w_i \in L^{\infty}(\mathbb{R}^N)$ such that $\displaystyle w_i := ess \inf_{x\in \mathbb{R}^N} w_i(x) > 0$. $M_i$ are bounded continuous   functions, $F$ belongs to $C^1(\mathbb{R}^{N}\times \mathbb{R}^{n})$ and satisfies adequate growth assumptions and   $F_{u_i}$ denotes the partial derivative of $F$ with respect to $u_i$ for all $1\leq i \leq n$. The functions $p_i(x), q_i(x)$ and $\gamma_i(x)$ are  continuous real-valued functions such that
	\begin{gather}
		1 <p_i^{-}\leq p_i(x) \leq p_i^{+}< q_i^{-}\leq q_i(x) \leq q_i^{+}
		< N,
	\end{gather}
	and
	$$\gamma_i(x) =(1-\mathcal{H}(k^3_{i})) p_i(x)+ \mathcal{H}(k^3_{i})q_i(x),$$

	for all $x\in \mathbb{R}^N$ with $k^3_{i}$ is given in $(\textbf{\textit{H}}_2)$, where $\displaystyle p_i^-:= \inf_{x\in \mathbb{R}^N}p_i(x)$,  $\displaystyle p_i^{+}:= \sup_{x\in \mathbb{R}^N}p_i(x)$, and analogously to $q_i^-, q_i^+, \gamma_i^-$ and $\gamma_i^+$, with $\displaystyle \mathcal{H}:\mathbb{R}_{0}^{+}\to \left\lbrace 0,1\right\rbrace  $ is given by
	\[
	\mathcal{H}(k_i)=\begin{cases}
		1 & \text{ if }        	k_i \geq 0,\\
		0 & \text{ if }            k_i<0.
	\end{cases}
	\]

	The operator $\mathcal{A}_{ji} : X_i\to \mathbb{R}^n$ with $j=1$ or $2$, and the operator  $ \mathcal{B}_i : X_i\to \mathbb{R}$, are respectively defined by
	\begin{equation}\label{AB}
		\mathcal{A}_{ji}(u_i)= a_{ji}(| u_i|^{p_i(x)}) | u_i|^{p_i(x)-2} u_i, \text{ and }\mathcal{B}_i(u_i)= \displaystyle\int_{\mathbb{R}^N }\dfrac{1}{p_i(x)}\Big(A_{1i}(|\nabla u_i|^{p_i(x)})+ w_i(x)A_{2i}(|u_i|^{p_i(x)})\Big)dx,
	\end{equation}
	where $X_i$ is the Banach space $X_i:= W_{w_i}^{1,p_i(x)}(\mathbb{R}^N)\cap W_{w_i}^{1,\gamma_i(x)}(\mathbb{R}^N)$, $A_{ji}(.)$ is the function $A_{ji}(t)=\displaystyle\int_{0}^{t}a_{ji}(k) dk$, and the function $a_i(.)$ is described in the hypothesis $(\textbf{\textit{H}}_1)$. \\

	In this article, we consider the function $ a_{ji}: \mathbb{R}^+ \to \mathbb{R}^+$  satisfying the following hypotheses for all $1\leq i \leq n$ and $j\in \{1,2\}$:

	\begin{itemize}
		\item[$(\textbf{\textit{H}}_1)$] The function $a_{ji}(.)$ is of   class $C^1$.

		\item[$(\textbf{\textit{H}}_2)$] There exist positive constants $k_{ji}^0, k_{ji}^1, k_{ji}^2$ and $k_{i}^3$ for all $1\leq i\leq n$, such that

		$$k_{ji}^0 + \mathcal{H}(k_{i}^3)k_{ji}^2 \tau^{\frac{q_i(x)-p_i(x)}{p_i(x)}} \leq a_{ji}(\tau ) \leq k_{ji}^1 + k_{i}^3\tau^{\frac{q_i(x)-p_i(x)}{p_i(x)}},
		$$
		for all $\tau \geq 0$ and for almost every $x\in \mathbb{R}^N.$

		\item[$(\textbf{\textit{H}}_3)$] There exist $c_{ji}>0$ such that

		$$ \min \left\lbrace a_{ji}(\tau^{p_i(x)})\tau^{p_i(x)-2}, a_{ji}(\tau^{p_i(x)})\tau^{p_i(x)-2}
		+\tau \frac{\partial(a_{ji}(\tau^{p_i(x)})\tau^{p_i(x)-2})}{\partial \tau } \right\rbrace \geq c_{ji}\tau^{p_i(x)-2},
		$$
		for almost every $x\in \mathbb{R}^N$ and for all $\tau>0$.

		\item[$(\textbf{\textit{H}}_4)$]
		There exists positive constants $\beta_{ji}$ for all $i\in \left\lbrace 1,2,\dots,n\right\rbrace $ such that

		$$ A_{ji}(\tau )\geq \frac{1}{\beta_{ji}}a_{ji}(\tau)\tau,$$
		for all $\tau \geq 0$.

		\item[$(\textbf{\textit{M}})$] $M_i: \mathbb{R}^+ \longrightarrow \mathbb{R}$ are continuous and increasing functions such that $0 < m_0=\displaystyle\min_{1 \leq i \leq n} M_i(0)  \leq M_i(t) \leq m_1=\max_{1 \leq i \leq n}M_i(t)$, for all $t \geq 0$, $i\in \left\lbrace 1,2,\dots,n\right\rbrace $. \\
	\end{itemize}

     Recently, a great attention has been focused on the study of this type of problems. They appear in mathematical models in different branches in science as electrorheological fluids \cite{Ru1,Ru2}, elastic mechanics \cite{Zhik1}, stationary thermorheological viscous flows of non-Newtonian fluids \cite{Raj}, image processing \cite{Chen1}, and mathematical description of the processes filtration of barotropic gas through a porous medium \cite{Anto2}. \\

     The system (\ref{s1.1}) is related (in the case of a single equation)  to a model firstly proposed by Kirchhoff in 1883 as the stationary version of the Kirchhoff equation

     \begin{equation}
     	\label{e1.1}
     	\rho\frac{\partial^2u}{\partial t^2}-\left( \dfrac{\rho_0}{h}+\dfrac{E}{2L}\int_{0}^{L}\left\vert\dfrac{\partial u(x)}{\partial x}\right\vert^2dx \right)\dfrac{\partial^2u}{\partial x^2}=0,
     \end{equation}
     where $\rho$, $\rho_0$, $E$ and $L$ are constants. This equation   extends the classical D'Alembert's wave equation by considering the effects of the changes in the length of the strings during the vibrations. A distinguishing feature of equation (\ref{e1.1}) is that the equation contains a nonlocal coefficient $\dfrac{\rho_0}{h}+\dfrac{E}{2L}{\displaystyle\int_{0}^{L}\left\vert\dfrac{\partial u}{\partial x}\right\vert^2dx}$ which depends on the average $\dfrac{1}{L}{\displaystyle\int_{0}^{L}\left\vert\dfrac{\partial u}{\partial x}\right\vert^2dx}$, and hence the equation is no longer a pointwise equation. The parameters in equation (\ref{e1.1}) have the following meanings: $E$ is the Young modulus of the material, $\rho$ is the mass density, $L$ is the length of the string, $h$ is the area of cross-section, and $\rho_0$ is the initial tension (see \cite{Kir}). \\

	Now, in order to illustrate the degree of generality of the kind of problems studied here, with adequate hypotheses on the functions $a_{ji}$, in the following we present more some examples of problems which are also interesting from the mathematical point of view and have a wide range of applications in physics and related sciences. \\

	\textbf{Example I.} Considering $a_{ji}\equiv 1$, we have that $a_{ji}$ satisfies the  $(\textbf{\textit{H}}_1),(\textbf{\textit{H}}_2)$ and $(\textbf{\textit{H}}_3)$ with $k^0_{ji}=k^1_{ji}=1$  and $k^2_{ji}>0$ and $k^3_{i}=0$ for all $i \in \{1,2,...,N\}$ and for $j=1$ or $2$. In this case we are studying problem:

	\begin{eqnarray}
		\label{ex1.1}
		\begin{cases}
			-M_i\Big(\mathcal{B}_i(u_i)\Big)\Big(\Delta_{p_i(x)}u_i -w_i(x)| u_i|^{p_i(x)-2} u_i\Big)
			= \lambda F_{u_i}(x,u)  & \text{ in } \mathbb{R}^N, \\
			\quad	u_i \in W_{w_i}^{1,p_i(x)}(\mathbb{R}^N)\\

		\end{cases}
	\end{eqnarray}

	where $\mathcal{B}_i(u_i)=\displaystyle\int_{\mathbb{R}^N }\dfrac{1}{p_i(x)}\left(|\nabla u_i|^{p_i(x)}+ w_i(x)|u_i|^{p_i(x)} \right)dx$. The operator $ \Delta_{p_i(x)}u_i:=\textrm{div}\,(|\nabla u_i|^{p_i(x)-2}\nabla u_i)$ is so-called $p_i(x)$-Laplacian, which coincides with the usual $p_i$-Laplacian when $p_i(x)=p_i$, and with the Laplacian when $p_i(x)=2$. \\

	\textbf{Example II.} Considering $a_{ji}(t)= 1+ t^{\frac{q_i(x)-p_i(x)}{p_i(x)}}$,
	we have that $a_{ji}$ satisfies the  $(\textbf{\textit{H}}_1),(\textbf{\textit{H}}_2)$ and $(\textbf{\textit{H}}_3)$ with $k^0_{ji}=k^1_{ji}=k^2_{ji}=\kappa_i^3=1$ for all $i\in \{1,2,...,n\}$ and for $j=1$ or $2$). In this case we are studying the $p_i\&q_i$-Laplacian equation:
	\begin{eqnarray}
		\label{ex1.12}
		\begin{cases}
			-M_i\Big(\mathcal{B}_i(u_i)\Big)\Big(\Delta_{p_i(x)}u_i +\Delta_{q_i(x)}u_i -w_i(x) (| u_i|^{p_i(x)-2} u_i +| u_i|^{q_i(x)-2} u_i) \Big)
			= \lambda F_{u_i}(x,u)  & \text{ in } \mathbb{R}^N, \\
			\qquad	u_i \in W_{w_i}^{1,p_i(x)}(\mathbb{R}^N)\cap W_{w_i}^{1,q_i(x)}(\mathbb{R}^N)\\
		\end{cases}
	\end{eqnarray}

	where $\mathcal{B}_i(u_i)=\displaystyle\int_{\mathbb{R}^N }\Bigg(\dfrac{1}{p_i(x)}\left(|\nabla u_i|^{p_i(x)}+ w_i(x)|u_i|^{p_i(x)} \right)+\dfrac{1}{q_i(x)}\left(|\nabla u_i|^{q_i(x)}+ w_i(x)|u_i|^{q_i(x)} \right)\Bigg)dx$.

	This class of problems comes, for example, from a general reaction-diffusion system

	\begin{equation} \label{DCE}
		u_t= \textrm{div}[D(u)\nabla u ]+ h(x,u),
	\end{equation}

	where $D(u)=|\nabla u|^{p(x)-2}+|\nabla u|^{q(x)-2}$, and the reaction term $h(x, u)$ is a polynomial of u with variable coefficients. This system has a wide range of applications in physics and related sciences, such as biophysics, plasma physics and chemical reaction design. In such applications, the function $u$ describes a concentration, the first term on the right-hand side of (\ref{DCE}) corresponds to the diffusion with a diffusion coefficient $D(u)$; whereas the second one is the reaction and relates to source and loss processes. Typically, in chemical and biological applications, (for further details, see \cite{Mahshid,He} references therein). \\

	We continued with other examples that are also interesting from mathematical point of view: \\

	\textbf{Example III.} Considering $a_{1i}(t)= 1+ \frac{t}{\sqrt{1+t^2}}$ and  $a_{2i}\equiv 1$, we have that $a_{1i}$ and $a_{2i}$ satisfies the  $(\textbf{\textit{H}}_1),(\textbf{\textit{H}}_2)$ and $(\textbf{\textit{H}}_3)$ with $k^0_{1i}=k^0_{2i}=k^1_{2i}=1$,$k^1_{1i}=2$, and $k^3_{i}=0$, $k^2_{1i}>0$ and $k^2_{2i}>0$. In this case we are studying problem:

	\begin{eqnarray}
		\label{ex1.13}
		\begin{cases}
			-M_i\Big(\mathcal{B}_i(u_i)\Big)\Big(\textrm{div}\,\Big( \Big(1+ \frac{|\nabla u_i|^{p_i(x)}}{\sqrt{1+|\nabla u_i|^{2p_i(x)}}}\Big)|\nabla u_i|^{p_i(x)-2}\nabla u_i\Big) -w_i(x)| u_i|^{p_i(x)-2} u_i\Big)
			= \lambda F_{u_i}(x,u)   & \text{ in } \mathbb{R}^N, \\
			\quad	u_i \in W_{w_i}^{1,p_i(x)}(\mathbb{R}^N),
		\end{cases}
	\end{eqnarray}

	where $\mathcal{B}_i(u_i)=\displaystyle\int_{\mathbb{R}^N }\frac{1}{p_i(x)}\Big(|\nabla u_i|^{p_i(x)} + \sqrt{1+|\nabla u_i|^{2p_i(x)}}+ w_i(x)|u_i|^{p_i(x)}\Big)dx$. \\

	The operator  $\textrm{div}\,\Big( \Big(1+ \frac{|\nabla u|^{p(x)}}{\sqrt{1+|\nabla u|^{2p(x)}}}\Big)|\nabla u|^{p(x)-2}\nabla u\Big)$ is so-called  $p(x)$-Laplacian like or so-called the generalized Capillary operator. The Capillarity can be briefly explained by considering the effects of two opposing forces: adhesion, i.e. the attractive (or repulsive) force between the molecules of the liquid and those of the container; and cohesion, i.e. the attractive force between the molecules of the liquid. The study of capillary phenomenon has gained much attention. This increasing interest is motivated not only by the fascination in naturally-occurring phenomena such as motion of drops, bubbles and waves, but also its importance in applied fields ranging from industrial and biomedical and pharmaceutical to microfluidic systems see \cite{Ni}. \\

	\textbf{Example IV.} Considering $\displaystyle a_{1i}(t)= 1+ t^{\frac{q_i(x)-p_i(x)}{p_i(x)}} +\frac{1}{(1+t)^{\frac{p_i(x)-2}{p_i(x)}}}$ and $\displaystyle a_{2i}(t)= 1+ t^{\frac{q_i(x)-p_i(x)}{p_i(x)}}$, we have that $\displaystyle a_{1i}$ and $\displaystyle a_{2i}$  satisfies the  $\displaystyle (\textbf{\textit{H}}_1),(\textbf{\textit{H}}_2)$ and $(\textbf{\textit{H}}_3)$ with $k^0_{1i}=k^0_{2i}=k^1_{2i}=1$, $k^1_{1i}= 2$ and $k^3_{i}=k^2_{1i}=k^2_{2i}=1$. In this case we are studying problem

	\begingroup\makeatletter\def \f@size{8}\check@mathfonts
	\begin{eqnarray}
		\label{ex1.13}
		\begin{cases}
			-M_i\Big(\mathcal{B}_i(u_i)\Big)\Bigg[\Delta_{p_i(x)}u_i +\Delta_{q_i(x)}u_i +\textrm{div}\,\Big( \dfrac{|\nabla u_i|^{p_i(x)-2}\nabla u_i}{(1+ |\nabla u_i|^{p_i(x)})^{\frac{p_i(x)-2}{p_i(x)}}}\Big)-w_i(x) \Big(| u_i|^{p_i(x)-2} u_i +| u_i|^{q_i(x)-2} u_i \Big)\Bigg]
			= \lambda F_{u_i}(x,u)  & \text{ in } \mathbb{R}^N, \\
			\qquad	u_i \in W_{w_i}^{1,p_i(x)}(\mathbb{R}^N)\cap W_{w_i}^{1,q_i(x)}(\mathbb{R}^N),
		\end{cases}
	\end{eqnarray}
	\endgroup

	where $\mathcal{B}_i(u_i)=\displaystyle\int_{\mathbb{R}^N }\Bigg(\dfrac{1}{p_i(x)}\left(|\nabla u_i|^{p_i(x)}+ w_i(x)|u_i|^{p_i(x)} \right)+\dfrac{1}{q_i(x)}\left(|\nabla u_i|^{q_i(x)}+ w_i(x)|u_i|^{q_i(x)} \right) + \frac{1}{2}(1+ |\nabla u_i|^{p_i(x)})^{\frac{2}{p_i(x)}}\Bigg)dx$. \\

	In the literature, the existence and multiplicity of the solutions for quasilinear elliptic systems have been studied by many authors (\cite{Afrouzi,BonaCand1,Bona1,Bona2,Chems0,Dai1,Dai2,Liu1, Shi1}), where the general variable exponent elliptic operators $\mathcal{A}_{ji}$ in divergence form and the nonlinear potential $F$ have different mixed growth conditions. For example in \cite{Djellit} the authors show the existence of nontrivial solutions for the following $(p,q)$-Laplacian system
	\begin{eqnarray}
		\label{ss1.2}
		\begin{cases}
			\displaystyle -\Delta _{p}u=\frac{\partial F}{\partial u}(x,u,v)\quad \text{in }\mathbb{R}^N \\
			& \\
			\displaystyle -\Delta _{q}v=\frac{\partial F}{\partial v}(x,u,v)\quad \text{in }\mathbb{R}^N
		\end{cases}
	\end{eqnarray}
	where the potential function $F$ satisfies mixed and subcritical growth conditions and, in addition, is supposed to be intimately connected with the first eigenvalue of the $(\Delta _{p},\Delta _{q})$-operator. They apply the Mountain Pass theorem to get existence of nontrivial solutions.  In \cite{Chems}, the authors using an abstract
	critical point result of Bonanno and Marano established the existence of an interval $\Lambda \subseteq [0, +\infty[$ such that for each $\lambda \in \Lambda$ the quasilinear elliptic system
	\begin{equation}
		\label{sSS1.1}
		-\Delta _{p_i(x)}u_i +a_i(x)|u_i|^{p_i(x)-2}u_i=\lambda F_{u_i}(x,u_1,u_2,...,u_n) \quad\text{in }\mathbb{R}^N,
	\end{equation}
	for $1\leq i \leq n$, where $a_i \in L^{\infty}(\mathbb{R}^N)$ such that $\displaystyle a_i := ess \inf_{x\in \mathbb{R}^N} a_i(x) > 0$ and  $\lambda $ is a positive parameter, has at least three distinct nontrivial solutions, and in \cite{Chems1} some similar results were obtained for the following classe of nonlocal elliptic systems
	 	\begin{eqnarray}
		\label{ssss1.1}
		\begin{cases}
			-M_1\left({\displaystyle\int_{\mathbb{R}^N }\dfrac{1}{p(x)}\left(|\nabla u|^{p(x)}+ a(x)|u|^{p(x)}\right)dx}\right)
			\big(\Delta _{p(x)}u -a(x)|u|^{p(x)-2}u\big)=\lambda F_{u}(x,u,v) &\quad\text{in }\mathbb{R}^N, \\

			-M_2\left({\displaystyle\int_{\mathbb{R}^N }\dfrac{1}{q(x)}\left(|\nabla v|^{q(x)}+ a(x)|v|^{q(x)}\right)dx}\right)
			\big(\Delta _{q(x)}v -b(x)|v|^{q(x)-2}v\big)=\lambda F_{v}(x,u,v) &\quad\text{in }\mathbb{R}^N;
		\end{cases}
	\end{eqnarray}

	where $a $, $b \in L^{\infty}(\mathbb{R}^N)$ such that $\displaystyle a := ess \inf_{x\in \mathbb{R}^N} a(x) > 0$ and $\displaystyle b := ess \inf_{x\in \mathbb{R}^N} b(x) > 0$. $M_1$ and $M_2$ are bounded continuous   functions,   $F$ belongs to $C^1(\mathbb{R}^{N}\times \mathbb{R}^{2})$ and verifies some mixed growth conditions.\\

    This paper's primary purpose is to establish the existence of some interval which includes $\lambda$, where the system (\ref{s1.1}) admits at least three weak solutions, which extend, complement and complete in several ways some of many works in particular the results in \cite{Chems,Chems1}, by means of an abstract critical points result of G. Bonanno  and S.A. Marano \cite{BonaMara}, which is a more precise version of Theorem 3.2 of \cite{BonaCand1}. For other basic notations and definitions we refer to \cite{Zeidler}. \\

	Now, and for the convenience of the readers, we recall the three critical points theorem of G. Bonanno and S.A. Marano \cite{BonaMara} which is our main tool to prove the results. Here, $X^*$ denotes the dual space of $X$.

	\begin{lemma}[see {\cite[Theorem 3.6]{BonaMara}}]\label{lem1.1}
		Let $X$ be a reflexive real Banach space; $\Phi:X\rightarrow \mathbb{R}$ be a
		coercive, continuously G\^{a}teaux differentiable and sequentially weakly
		lower semicontinuous functional whose G\^{a}teaux derivative admits a
		continuous inverse on $X^*$; $\Psi:X\rightarrow \mathbb{R}$ be a continuously G\^{a}teaux differentiable functional whose G\^{a}teaux derivative is
		compact such that
		$$
		\Phi(0)=\Psi(0)=0.
		$$
		Assume that there exist $r>0$ and $\overline{x}\in X,$ with $r<\Phi
		(\overline{x}),$ such that
		\begin{enumerate}
			\item[$\rm(a_1)$] $\displaystyle \frac{\displaystyle \sup_{\Phi(\bar x)\leq r}\Psi(\bar x)}{r}<\frac{\Psi
				(\overline{x})}{\Phi(\overline{x})}$;
			\item[$\rm(a_2)$] for each $\displaystyle \lambda\in \Lambda_r:=]\frac{\Phi (\overline{x})}{\Psi(\overline{x})},\frac{r}{\displaystyle \sup_{\Phi(\bar x)\leq r}\Psi(\bar x)}[$ the
			functional $\Phi-\lambda\Psi$ is coercive.
		\end{enumerate}
		Then, for each $\lambda\in\Lambda_r$ the functional $\Phi-\lambda
		\Psi$ has at least three distinct critical points in $X$.
	\end{lemma}

	\textit{Organization of the paper.} The rest of the paper is organized as follows: Section 2 contains some basic preliminary knowledge of the variable exponent spaces and some results which will be needed later. Finally, in Section 3, we state and establish our main result.

	\section{Preliminaries and basic notations}\label{sec2}
	In this section, we introduce some definitions and results which will be used in the next section. Firstly, we introduce some theories of Lebesgue-Sobolev spaces with variable exponent. The details can be found in \cite{Dien1, Fan3, Kov1, Rad}. Denote by $S(\mathbb{R}^N)$ the set of all measurable real functions on $\mathbb{R}^N$. Set

	\begin{equation*}
		C_{+}(\mathbb{R}^N)= \{  p \in C(\mathbb{R}^N) :  \inf_{x\in \mathbb{R}^N}p(x)> 1 \} .
	\end{equation*}

	For any $ p \in C_{+}(\mathbb{R}^N)$ we define
	\begin{gather*}
		p^{-}:= \inf_{x\in \mathbb{R}^N}p(x) \quad\text{ and} ~p^{+}:= \sup_{x\in \mathbb{R}^N}p(x) .
	\end{gather*}

	For any $p \in C_{+}(\mathbb{R}^N)$, we define the variable exponent Lebesgue space as

	\begin{gather*}
		L^{p(x)}(\mathbb{R}^N)= \left\{ u \in S(\mathbb{R}^N) : \int_{\mathbb{R}^N}|u(x)|^{p(x)}dx <\infty \right\},
	\end{gather*}

	endowed with the Luxemburg norm

	\begin{equation*}
		|u|_{p(x)}:=|u|_{L^{p(x)}(\mathbb{R}^N)} = \inf \left\{ \mu >0 : \int_{\mathbb{R}^N }\left|\frac{u(x)}{\mu}\right|^{p(x)} dx \leq 1 \right \}.
	\end{equation*}

	Let $w \in S(\mathbb{R}^N)$, and $w(x)>0$ for a.e $x\in \mathbb{R}^N$. Define the weighted variable exponent Lebesgue space $L_{w}^{p(x)}(\mathbb{R}^N)$ by

	\begin{gather*}
		L_{w}^{p(x)}(\mathbb{R}^N)= \left\{ u \in S(\mathbb{R}^N) : \int_{\mathbb{R}^N}w(x)|u(x)|^{p(x)}dx <\infty \right\},
	\end{gather*}

	with the norm

	\begin{equation*}
		|u|_{p(x),w(x)}:=|u|_{L_w^{p(x)}(\mathbb{R}^N)} = \inf \left\{ \mu >0 : \int_{\mathbb{R}^N }w(x)\left|\frac{u(x)}{\mu}\right|^{p(x)} dx \leq 1 \right \}.
	\end{equation*}

	From now on, we suppose that $w \in L^{\infty}(\mathbb{R}^N)$ with $\displaystyle w := ess \inf_{x\in \mathbb{R}^N} w(x) > 0$. Then obviously $L_w^{p(x)}(\Omega)$ is a Banach space (see \cite{Cruz1} for details). \\

	On the other hand, the variable exponent Sobolev space $W^{1,p(x)}(\mathbb{R}^N)$ is defined by

	\[
	W^{1,p(x)}(\mathbb{R}^N)=\{ u\in L^{p(x)}(\mathbb{R}^N):| \nabla
	u| \in L^{p(x)}(\mathbb{R}^N)\},
	\]

	and is endowed with the norm

	\[
	\| u\| _{1,p(x)}:=\| u\|_{W^{1,p(x)}(\mathbb{R}^N)}
	=|u| _{p(x)}+| \nabla u| _{p(x)},   ~~~~ \forall u\in W^{1,p(x)}(\mathbb{R}^N ).
	\]

	Next, the weighted-variable exponent Sobolev space $W_{w}^{1,p(x)}(\mathbb{R}^N)$ is defined by

	\[
	W_{w}^{1,p(x)}(\mathbb{R}^N)=\{ u\in L_{w}^{p(x)}(\mathbb{R}^N):| \nabla
	u| \in L_{w}^{p(x)}(\mathbb{R}^N)\},
	\]

	with the norm

	\begin{equation*}
		\| u\| _{p(x),w}:= \inf \left\{ \mu >0 :  \int_{\mathbb{R}^N } \left|\frac{\nabla u(x)}{\mu}\right|^{p(x)} +w(x)\left|\frac{u(x)}{\mu}\right|^{p(x)} dx \leq 1 \right \},  \forall u \in W_{w}^{1,p(x)}(\mathbb{R}^N).
	\end{equation*}

	Then the norms $\| u\| _{p(x),w}$ and $\| u\| _{p(x)}$ are equivalent in $W_{w}^{1,p(x)}(\mathbb{R}^N)$.

	If $p^->1$, then the spaces $L^{p(x)}(\mathbb{R}^N)$, $W^{1,p(x)}(\mathbb{R}^N)$ and $W_w^{1,p(x)}(\mathbb{R}^N)$ are separable, reflexive and uniformly convex Banach spaces. \\

	Here we display some facts which will be used later.

	\begin{proposition}[see \cite{Dien1, Fan3}] \label{prop1}
			The topological dual space of
		$L^{p(x)}(\Omega)$ is $L^{p'(x)}(\Omega)$, where
		\[
		\frac{1}{p(x)}+\frac{1}{p'(x)}=1.
		\]
		Moreover, for any $(u,v)\in L^{p(x)}(\Omega)\times L^{p'(x)}(\Omega)$,
		we have
		\[
		\Big| \int_{\Omega}uvdx\Big|
		\leq (\frac{1}{p^{-}}+ \frac{1}{(p')^{-}})| u|_{p(x)}| v| _{p'(x)}
		\leq 2| u| _{p(x)}| v|_{p'(x)}.
		\]
	\end{proposition}

	\begin{proposition}[see \cite{Dien1, Fan3}] \label{prop2}
		Denote $\displaystyle \rho_{p(x)} (u):=\int_{\mathbb{R}^N}| u| ^{p(x)}dx$,  for all $u\in L^{p(x)}(\mathbb{R}^N)$. We have
		\[
		\min \{ | u| _{p(x)}^{p^{-}},| u| _{p(x)}^{p^{+}}\}
		\leq \rho_{p(x)} (u)\leq \max \{ | u| _{p(x)}^{p^{-}},| u| _{p(x)
		}^{p^{+}}\},
		\]

		and the following implications are true
		\begin{itemize}
			\item[(i)]  $|u| _{p(x)}<1$ (resp. $=1, >1$) $\Leftrightarrow \rho_{p(x)} (u)<1$
			(resp. $=1,>1$),

			\item[(ii)] $|u| _{p(x)}>1 \Rightarrow | u| _{p(x)}^{p^{-}}\leq \rho_{p(x)} (u)
			\leq | u| _{p(x)}^{p^{+}}$,

			\item[(iii)] $|u|_{p(x) }<1\Rightarrow | u| _{p(x)}^{p^{+}}\leq \rho_{p(x)} (u)
			\leq | u| _{p(x)}^{p^{-}}$.
		\end{itemize}
	\end{proposition}

Denote $\displaystyle \rho_{p(x),w}(u):=\int_{\mathbb{R}^N }\left( \left|\nabla u(x)\right|^{p(x)} +w(x)\left|u(x)\right|^{p(x)} \right)dx$,  for all $u \in W_{w}^{1,p(x)}(\mathbb{R}^N)$. From Proposition \eqref{prop2}, we have

		\begin{equation}
			\label{inq1}
			\| u\| _{p(x),w}^{p^{-}}  \leq \rho_{p(x),w}(u) \leq  \| u\| _{p(x),w}^{p^{+}}    \quad \text{ if} ~~~~ \| u\| _{p(x),w} \geq 1,
		\end{equation}
		\begin{equation}
			\label{inq2}
			\| u\| _{p(x),w}^{p^{+}}  \leq \rho_{p(x),w}(u) \leq  \| u\| _{p(x),w}^{p^{-}}    \quad \text{ if} ~~~~ \| u\| _{p(x),w} \leq 1.
		\end{equation}

	\begin{proposition}[see \cite{Edmu1}] \label{prop3}
		Let $p(x)$ and $q(x)$ be measurable functions such that $p\in L^{\infty }(\mathbb{R}^N)$ and $1\leq p(x), q(x)< \infty $ almost everywhere in
		$\mathbb{R}^N$. If $u\in L^{q(x)}(\mathbb{R}^N)$, $u\neq 0$, then we have

		\begin{gather*}
			| u| _{p(x)q(x)}\leq 1\Rightarrow |u| _{p(x)q(x)}^{p^{-}}
			\leq \big|| u| ^{p(x)}\big| _{q(x)}\leq | u| _{p(x)q(x)}^{p^{+}},
			\\
			| u| _{p(x)q(x)}\geq 1\Rightarrow |u| _{p(x)q(x)}^{p^{+}}
			\leq \big| | u| ^{p(x)}\big| _{q(x)}\leq | u| _{p(x)q(x)}^{p^{-}}.
		\end{gather*}

		In particular, if $p(x)=p$ is constant, then

		\[
		| | u| ^{p}| _{q(x)}
		=|u| _{pq(x)}^{p}.
		\]
	\end{proposition}

	Now, for all $x\in \mathbb{R}^N$, denote
	\[
	p^{\partial }(x)=\begin{cases}
		\frac{N p(x)}{N-p(x)} &\text{for } p(x)<N \\
		+\infty &\text{for }p(x)\geq N
	\end{cases}
	\]
	the critical Sobolev exponent of $p(x)$.  We have

	\begin{proposition}[see \cite{Dien1,Edmu1}] \label{prop4}
		Let $r\in C_{+}^{0,1}(\mathbb{R}^N)$, the space of Lipschitz-continuous
		functions defined on $\mathbb{R}^N$. Then, there exists a positive constant $c_{r}$ depending on $r$ such that

		\[
		| u| _{r^{\partial }(x)}\leq c_{r}\|u\|_{r(x),w},   \quad \forall
		u\in W_w^{1,r(x)}(\mathbb{R}^N).
		\]
	\end{proposition}

	\begin{proposition}[see \cite{Dien1,Edmu1}] \label{prop5}
		Assume that $p\in  C(\mathbb{R}^N) $ with $p(x)>1$ for each $x\in \mathbb{R}^N$. If $r\in C(\mathbb{R}^N)$ is such that $1<r(x)<p^{\partial}(x)$ for each $x\in \Omega$, then there exists a continuous and compact embedding $W^{1,p(x)}(\mathbb{R}^N)\hookrightarrow L^{r(x)}(\mathbb{R}^N )$.
	\end{proposition}

	In the following, we shall use the product space

	\[
	X:=\prod_{i=1}^{n}\Big(W_{w_i}^{1,p_i(x)}(\mathbb{R}^N )\cap W_{w_i}^{1,\gamma_i(x)}(\mathbb{R}^N )\Big),
	\]

	equipped with the norm

	\[
	\| u\|:=\max \left\lbrace  \|  u_i\|_{w_i} \right\rbrace ,~~~~~ \forall u=(u_1, \dots, u_n)\in X,
	\]

	where $\| u_i\|_{w_i}:= \|  u_i\|_{p_i(x),w_i}+ \mathcal{H}(\kappa_i^3)\| u_i\|_{q_i(x),w_i}$  is the norm in $W_{w_i}^{1,p_i(x)}(\mathbb{R}^N )\cap W_{w_i}^{1,\gamma_i(x)}(\mathbb{R}^N )$. We denote  $X^{\star}$ the topological dual of $X$ equipped with the usual dual norm.\par

	\begin{definition}
		Let $X$ be a Banach space, an element $u=(u_1,u_2,...,u_n) \in X$ is called a weak solution of the system \eqref{s1.1} if

		\begin{center}
			$\displaystyle \sum_{i=1}^{n} M_i\left(\mathcal{B}_i(u_i)\right)\int_{\mathbb{R}^N} \Big(\mathcal{A}_{1i}(\nabla u_i) \nabla v_i + w_i(x)\mathcal{A}_{2i}(u_i) v_i \Big)\,dx - \sum_{i=1}^{n}\displaystyle\int_{\Omega} \lambda F_{u_i}(x,u_1,...u_n)v_i\,dx =0,$
		\end{center}

		for all $\displaystyle v=(v_1,v_2,...,v_n)\in X= \prod_{i=1}^{n}(W_{w_i}^{1,p_i(x)}(\mathbb{R}^N )\cap W_{w_i}^{1,\gamma_i(x)}(\mathbb{R}^N ))$.
	\end{definition}

	For every $u = (u_1, \dots, u_n)$ in $X$, let us define the functional $\Psi$ by

	\begin{center}
		$\displaystyle \Psi(u) = \int_{\mathbb{R}^N} F(x, u_1, \dots, u_n) dx$.
	\end{center}

	Under the assumptions below $(\mathcal{F}_1)$ and $(\mathcal{F}_2)$, we have $\displaystyle \Psi \in C^1(X, \mathbb{R})$ and its derivative is given by

	\begin{center}
		$\displaystyle \Psi^\prime(u) v = \sum_{i = 1}^n D_i \Psi(u) v_i$,
	\end{center}

	where

	\begin{center}
		$\displaystyle D_i \Psi(u) v_i = \int_{\mathbb{R}^N} \frac{\partial F}{\partial u_i}(x, u_1, \dots, u_n) v_i dx$,
	\end{center}

for all $\displaystyle v=(v_1,v_2,...,v_n)\in X$. Define also  $\Phi$ in $X$  by

	\begin{center}
		$\displaystyle \Phi(u)=\sum_{i=1}^{n}\Phi_i(u_i)$,
	\end{center}
where
$$\Phi_i(u_i)=
\widehat{M_i}\left(\mathcal{B}_i(u_{i}(x))\right),$$
for any $u=(u_1,\dots,u_n)$ in $X$, with $\displaystyle \widehat{M}_i(t):=\int_{0}^{t}M_i(s)ds$ for all $t\geq0$, ($i\in \{1,\dots,n\}$).

	We recall that $\Phi$ is a $C^1$-functional, weakly lower semi-continuous and its derivative is given by

	\begin{center}
		$\displaystyle \Phi^\prime (u) v =  \sum_{i = 1}^n D_i \Phi(u) v_i$,
	\end{center}

	where

	\begin{center}
		$\displaystyle D_i \Phi(u) v_i = M_i\left(\mathcal{B}_i(u_i)\right)\int_{\mathbb{R}^N} \Big(\mathcal{A}_{1i}(\nabla u_i) \nabla v_i + w_i(x)\mathcal{A}_{2i}(u_i) v_i \Big)\,dx$,
	\end{center}

	for all $v = (v_1, \dots, v_n) \in X$. \\

	The Euler-Lagrange functional associated to the system \eqref{s1.1} is defined by

	\begin{center}
		$\displaystyle E_{\lambda}(u):= \Phi(u)-\lambda \Psi(u) \quad \forall u \in X$.
	\end{center}

	Since $E_{\lambda} \in C^{1}(X, \mathbb{R})$ and for all $v = (v_1, \dots, v_n) \in X$, we have

	\begin{center}
		$\displaystyle E_{\lambda}^\prime(u) v = \Phi^\prime (u) v - \lambda \Psi^\prime (u) v$
	\end{center}

	Consequently, $u \in X$ is a weak solution of \eqref{s1.1} if and only if $u$ is a critical point of $E_\lambda$. \\

	We end this section by stating some hypotheses that we will be using later. \\

		\subsection*{\textbf{Hypotheses}} We assume some growth conditions:
	\begin{itemize}
		\item[$(\mathcal{F}_1)$] $F\in C^{1}(\mathbb{R}^N\times \mathbb{R}^n,\mathbb{R})$ and $F(x,0,\dots, 0)=0$.

		\item[$(\mathcal{F}_2)$] There exist positive functions $b_{ij}$ ($1\leq i,j \leq n$), such that
		\begin{gather*}
			\Big| \frac{\partial F}{\partial u_i}(x,u_1,...,u_n)\Big| \leq
			\sum_{j=1}^{n} b_{ij}(x)| u_j|^{\mu_{ij}-1},
		\end{gather*}
		where
		$\displaystyle 1<\mu_{ij}<\inf_{x\in \mathbb{R}^N}\gamma_i(x)$ for all $x\in \mathbb{R}^N$ and for all $i\in \left\lbrace 1,2,...,n\right\rbrace$. The weight-functions $b_{ii}$ (resp  $b_{ij}$ if $i\neq j$) belong to the generalized Lebesgue spaces $L^{\alpha _i}(\mathbb{R}^N)$ (resp $L^{\alpha_{ij} }(\mathbb{R}^N)$), with
		\[
		\alpha_i(x)=\frac{\gamma_i(x)}{\gamma_i(x)-1}, \quad \alpha_{ij} (x)
		=\frac{\gamma_i^{\partial }(x)\gamma_j^{\partial }(x)}{\gamma_i^{\partial }(x)\gamma_j^{\partial }(x)-\gamma_i^{\partial
			}(x)-\gamma_j^{\partial }(x)}.
		\]

	\textbf{Example:}
	We give an  example of potential  $F$ satisfying hypotheses $(\mathcal{F}_1)$ and $(\mathcal{F}_2)$ for $n=2$. Let
	\[
	F(x,u_1,u_2)=a(x)|u_1|^{l_1(x)}|u_2|^{l_2(x)},
	\]
	where $\displaystyle \frac{l_1(x)}{\gamma_1(x)}+\frac{l_2(x)}{\gamma_2(x)}< 1$ and $a$ is a positive function in $L^{e(x)}(\mathbb{R}^N ) $ such that

	\begin{center}
	    $\displaystyle e(x)= \dfrac{\gamma_1^{\partial}(x)\gamma_2^{\partial}(x)}{\gamma_1^{\partial}(x)\gamma_2^{\partial}(x) -l_1(x) \gamma_2^{\partial}(x)-l_2(x) \gamma_2^{\partial}(x)}$,
	\end{center}

	for each $x\in \mathbb{R}^N $. \\

	We can easily verify that $F(x,u_1,u_2)$ satisfies the condition $(\mathcal{F}_1)$. Moreover, by using Young inequality we easily check that the condition $(\mathcal{F}_2)$. \\

		\item[$(\mathcal{F}_3)$] Assume that there exist $r>0$ and $z=(z_1,...,z_n)\in X$ such that the following conditions are satisfied:
		\begin{itemize}
			\item[$(C_1)$] $\displaystyle \sum_{i=1}^{n}\frac{m_0^{\star}}{p_i^{+}}\Big(\min \left\lbrace \| z_i\|_{p_i(x),w_i}^{p_i^-} ,\| z_i\|_{p_i(x),w_i}^{p_i^+} \right\rbrace + \mathcal{H}(k_i^3) \min \left\lbrace \| z_i\|_{q_i(x),w_i}^{q_i^-} ,\| z_i\|_{q_i(x),w_i}^{q_i^+} \right\rbrace\Big)>r.$

				\begingroup\makeatletter\def \f@size{9}\check@mathfonts
			\item[$(C_2)$] $\displaystyle \frac{\displaystyle \int_{\mathbb{R}^N} \sup_{(\xi_1,\dots,\xi_n) \in K(\frac{sr}{m_0^{\star}})}F(x,\xi_1,\dots,\xi_n)dx }{r} < \dfrac{\displaystyle \int_{\mathbb{R}^N} F(x,z_1, \dots ,z_n) dx}{\displaystyle m_1 \sum_{i=1}^n \Big(\max \left\lbrace \| z_i\|_{p_i(x),w_i}^{p_i^-} ,\| z_i\|_{p_i(x),w_i}^{p_i^+} \right\rbrace + \mathcal{H}(k_i^3) \max \left\lbrace \| z_i\|_{q_i(x),w_i}^{q_i^-} ,\| z_i\|_{q_i(x),w_i}^{q_i^+} \right\rbrace\Big)}$
				\endgroup
		\end{itemize}

		where

		\begin{equation}\label{defK}
			\displaystyle K(t):= \left\lbrace (\xi_1,...\xi_n)\in \mathbb{R}^N: \displaystyle \sum_{i=1}^{n} \bigg(\min \left\lbrace |\xi_i|_{p_i^{\partial}(x)}^{{(p_i^{\partial})^{-}}}, |\xi_i|_{p_i^{\partial}(x)}^{{(p_i^{\partial})^{+}}} \right\rbrace + \mathcal{H}(k_i^3) \min \left\lbrace |\xi_i|_{q_i^{\partial}(x)}^{{(q_i^{\partial})^{-}}}, |\xi_i|_{q_i^{\partial}(x)}^{{(q_i^{\partial})^{+}}} \right\rbrace \bigg)\leq t \right\rbrace,
		\end{equation}
		and
		\begingroup\makeatletter\def \f@size{9}\check@mathfonts
		\begin{equation}\label{defs}
				m_0^{\star} = \displaystyle m_0 \min_{1 \leq i \leq n}\Big\{\displaystyle \frac{\min\big\{\displaystyle \min\{k_{1i}^0,k_{2i}^0\}, \displaystyle \min \{k_{1i}^2,k_{2i}^2\}\big\}}{\displaystyle \max\{\beta_{1i}, \beta_{2i}\}}\Big\}
				\text{ and }
				\displaystyle s = \max \left\lbrace p_i^+ \min_{1\leq i\leq n} \left\lbrace c_{p_i(x)}^{{p_i^{-}}}, c_{p_i(x)}^{{p_i^{+}}}\right\rbrace , p_i^+ \min_{1\leq i\leq n} \left\lbrace c_{q_i(x)}^{{q_i^{-}}}, c_{q_i(x)}^{{q_i^{+}}}\right\rbrace \right\rbrace ,
		\end{equation}
	\endgroup
		with $t>0$, such that $c_{p_i(x)}$ and $c_{q_i(x)}$ representing the constants defined in Proposition \ref{prop4}.
	\end{itemize}

	\section{Main results}

	The following lemmas are needed in the proof of our main results, namely Theorem \ref{Th3.3} below.

	\begin{lemma}\label{lem.3.1}
		The functional $\Phi$ is continuously GÃ¢teaux differentiable and sequentially weakly lower semi-continuous, coercive and its GÃ¢teaux derivative admits a continuous inverse on $X^{\star}$.
	\end{lemma}

	\begin{proof}
		It is well known that the functional $\Phi$ is well defined and is a continuously GÃ¢teaux differentiable functional whose derivative at the point $u=(u_1,\dots,u_n) \in X$ is the functional $\Phi^\prime(u)$ given by
		\begin{center}
			$\displaystyle \langle \Phi^{\prime}\left(u\right),\varphi \rangle =\sum_{i=1}^{n}\langle \Phi_{i}^{\prime}\left(u_{i}\right), \varphi_{i} \rangle,$
		\end{center}

  where
				\begin{center}
			$\displaystyle \langle \Phi_{i}^\prime(u_i),\varphi_i \rangle = M_i\Big(\mathcal{B}_i(u_i)\Big) \int_{\mathbb R^N} \Big(\mathcal{A}_{1i}(\nabla u_i)\nabla \varphi_i + w_i(x)\mathcal{A}_{2i}(u_i)\varphi_i \Big) dx$,
		\end{center}
		for every $\varphi=(\varphi_1,...,\varphi_n) \in X$ and for all $i \in \{1,\dots,n\}$. Let us show that $\Phi$ is coercive. By using
		$(\textbf{H}_2)$, $(\textbf{H}_4)$ and $(\textbf{M})$, we have for all $u=(u_1,\dots,u_n) \in X$,

		\begin{center}
			$\displaystyle \begin{aligned} \Phi(u) & = \sum_{i=1}^n \Phi_{i}(u_i) = \sum_{i=1}^n  \widehat{M}_i\Big(\mathcal{B}_i(u_i)\Big) \\
				 & \geq \sum_{i=1}^n \frac{m_0}{p_i^+} \int_{\mathbb R^N} \Big(\mathcal{A}_{1i}(|\nabla u_i|) + w_i(x)\mathcal{A}_{2i}(|u_i|) \Big) dx \\
				 & \geq \sum_{i=1}^n \frac{m_0}{p_i^+} \int_{\mathbb R^N} \Big[\frac{1}{\beta_{1i}} a_{1i}(|\nabla u_i|^{p_i(x)})|\nabla u_i|^{p_i(x)} + \frac{1}{\beta_{2i}} w_i(x)a_{2i}(|u_i|^{p_i(x)})|u_i|^{p_i(x)} \Big] dx
				  \\ & \geq \sum_{i=1}^n  \frac{m_0}{p_i^+\max\{\beta_{1i}, \beta_{2i}\}} \int_{\mathbb R^N} \Big[a_{1i}(|\nabla u_i|^{p_i(x)})|\nabla u_i|^{p_i(x)} + w_i(x)a_{2i}(|u_i|^{p_i(x)})|u_i|^{p_i(x)} \Big] dx \\
				   & \geq \sum_{i=1}^n  \frac{m_0\min\{\min\{k_{1i}^0,k_{2i}^0\}, \min\{k_{1i}^2,k_{2i}^2\}\}}{p_i^+\max\{\beta_{1i}, \beta_{2i}\}} \cdot \Big(\rho_{p_i(x)w_i}(u_i) + \mathcal{H}(k_{i}^3) \cdot\rho_{q_i(x)w_i}(u_i) \Big),
				 \end{aligned}$
			\end{center}

		Hence, by inequalities \eqref{inq1} and $\eqref{inq2}$, we obtain
	\begin{center}
		$\displaystyle \begin{aligned} \Phi(u) & \geq \sum_{i=1}^n \frac{m_0\min\{\min\{k_{1i}^0,k_{2i}^0\}, \min\{k_{1i}^2,k_{2i}^2\}\}}{p_i^+\max\{\beta_{1i}, \beta_{2i}\}} \cdot \Big(\min\{\| u_i\|_{p_i(x),w_i}^{p_i^-} , \| u_i\|_{p_i(x),w_i}^{p_i^+}\}+ \mathcal{H}(k_{i}^3) \cdot \min\{\| u_i\|_{q_i(x),w_i}^{q_i^-}, \| u_i\|_{q_i(x),w_i}^{q_i^+}\} \Big). \end{aligned}$
		\end{center}

		This shows that $\Phi(u) \longrightarrow +\infty$ as $\| u\| \longrightarrow +\infty$ that is, $\Phi$ is coercive on $X$. Now, in order to show that the functional $\Phi^\prime: X \longrightarrow X^{\star}$ is strictly monotone, it suffices to prove that $\Phi$ is strictly convex. \newline

For all $i \in \{1,\dots,n\}$, the functional $\mathcal{B}_i : X_i \longrightarrow \mathbb{R}$ defined in \eqref{AB} is clearly a GÃ¢teaux differentiable at any $u_i \in X_i$, and its derivative is given by

		\begin{align*}
		 \displaystyle\langle \mathcal{B}_i^\prime(u_i),\varphi_i \rangle
			&= \int_{\mathbb R^N} \Big(\mathcal{A}_{1i}(|\nabla u_i|^{p_i(x)})\nabla \varphi_i + w_i(x)\mathcal{A}_{2i}(|u_i|^{p_i(x)})\varphi_i \Big) dx\\
			&= \int_{\mathbb R^N}\Big( a_{1i}(| \nabla u_i|^{p_i(x)}) | \nabla u_i|^{p_i(x)-2} \nabla u_i .\nabla\varphi_i
			+ w_i(x)a_{2i}(| u_i|^{p_i(x)}) | u_i|^{p_i(x)-2} u_i\varphi_i\Big)dx,
		\end{align*}
for all $\varphi_i\in X_i$. \\

	Taking into account the elementary inequalities (see, e.g., Auxiliary Results in \cite{Hurtado}) for any $ \varrho , \zeta \in \mathbb{R}^N $
	\begin{equation}\label{ineq1}
		\Big( a_{ji}(|\varrho|^{p_i(x)})|\varrho|^{p_i(x)-2}\varrho  - a_{ji}(|\zeta|^{p_i(x)})|\zeta|^{p_i(x)-2}\zeta
		\Big) \cdot\left(\varrho -\zeta\right)
		\geq \begin{cases}
			C_{p_i}|\varrho  -\zeta|^{p_i(x)} &\text{if } p_i(x)\geq2\\[4pt]
			C_{p_i}\frac{	|\varrho -\zeta|^{2}}{(|\varrho| +|\zeta| )^{p_i(x)-2}},\;(\varrho,\zeta)\neq(0,0) &\text{if }1<p_i(x)<2,
		\end{cases}
	\end{equation}

	where $\cdot$ denotes the standard inner product  in $\mathbb{R}^N$. Therefore, we have

		\begin{center}
			$\displaystyle \langle 	\mathcal{B}_i^\prime(u_i) - 	\mathcal{B}_i^\prime(v_i),u_i - v_i \rangle > 0$,
		\end{center}
for all  $u_i\neq v_i\in X_i$,
which means that $ \mathcal{B}_i'$ is strictly monotone. So, by  \cite[ Proposition 25.10]{Zeidler},
$\mathcal{B}_i$ is strictly convex. Moreover, since  the Kirchhoff function $M_i$ is nondecreasing, $\widehat{M}_i$ is convex in
$[0,+\infty[$. Thus, for every $u_i, v_i\in X_i$ with $u_i\neq v_i$, and every
$s,t\in (0,1)$ with $s+t=1$, we have
$$
\widehat{M}_i(\mathcal{B}_i(su_i+tv_i))<\widehat{M}_i(s\mathcal{B}_i(u_i)+t\mathcal{B}_i(v_i))\leq
s\widehat{M}_i(\mathcal{B}_i(u_i))+t\widehat{M}_i(\mathcal{B}_i(v_i)).
$$
	This shows that $\Phi_{i}$ is strictly convex in $ W_{w_i}^{1, p_i(x)}\left(\mathbb{R}^N\right) \cap W^{1, \gamma_i(x)}_{w_i}\left(\mathbb{R}^N\right)$ for all $i \in \{1,\dots,n\}$. Hence, $\Phi$ is strictly convex in $X$, and therefore $\Phi^\prime= \displaystyle\sum_{i=1}^{n} \Phi_{i}^\prime$ is strictly monotone.\\

	It's clear that $\Phi^\prime$ is an injection since $\Phi^\prime$ is strictly monotone operator in $X$.\newline
	 Moreover, since we have

		\begin{center}
			$\displaystyle \begin{aligned} \displaystyle \lim_{\|u\| \longrightarrow + \infty} \frac{\langle \Phi^\prime(u),u \rangle}{\|u\|} & = \displaystyle \lim_{\|(u_1,...,u_n)\| \to + \infty} \frac{\displaystyle\sum_{i=1}^{n} M_i\left(\mathcal{B}_i(u_i)\right)\int_{\mathbb{R}^N} \Big(\mathcal{A}_{1i}(\nabla u_i) \nabla u_i + w_i(x)\mathcal{A}_{2i}(u_i) u_i \Big)\,dx}{\|(u_1,...,u_n)\|}\\
				 & \geq \displaystyle \lim_{\|(u_1,...,u_n)\| \to + \infty} \frac{m_0 \displaystyle\sum_{i=1}^{n} \int_{\mathbb{R}^N} \Big(a_{1i}(| \nabla u_i|^{p_i(x)}) | \nabla u_i|^{p_i(x)} + w_i(x)a_{2i}(| u_i|^{p_i(x)}) |  u_i|^{p_i(x)} \Big)\,dx}{\|(u_1,...,u_n)\|}\\
				  & \geq \displaystyle \lim_{\|(u_1,...,u_n)\| \to+ \infty} \frac{m_0\displaystyle\min_{1\leq i\leq n}\Big\{\min\{k_{1i}^0,k_{2i}^0\}, \min\{k_{1i}^2,k_{2i}^2\}\Big\}\sum_{i=1}^{n} \Big( \| u_i\|_{p_i(x),w_i}^{p_i^+} + \mathcal{H}(k_{i}^3) \cdot \| u_i\|_{q_i(x),w_i}^{q_i^+} \Big) }{ \displaystyle\max_{1\leq i\leq n} \left\lbrace 	\|  u_i\|_{p_i(x),w_i}+ \mathcal{H}(\kappa_i^3)\| u_i\|_{q_i(x),w_i}\right\rbrace} =+\infty.
			 \end{aligned}$
		\end{center}

		Then, we deduce that $\Phi^\prime$ is coercive. Thus $\Phi^\prime$ is a surjection. Now, since $\Phi^\prime$ is semicontinuous in $X$, then by applying Minty-Browder theorem (Theorem 26.A of \cite{Zeidler}), we conclude that $\Phi^\prime$ admits a continuous inverse on $X^{\star}$. Moreover, the monotonicity of $\Phi^\prime$ on $X^{\star}$ ensures that $\Phi$ is sequentially lower semi-continuously on $X$ (see \cite{Zeidler}, Proposition 25. 20). The proof of the lemma is complete.
	\end{proof}

	\begin{lemma}\label{lem.3.2}
		Under assumptions $(\mathcal{F}_1)$ and $(\mathcal{F}_2)$, the functional $\Psi$ is well defined and is of class $C^1$ on $X$. Moreover, its derivative is given by

		\begin{center}
			$\displaystyle \displaystyle\langle\Psi^\prime(u), \varphi\displaystyle\rangle = \sum_{i=1}^n \frac{\partial F}{\partial u_i}(x,u) \varphi_i$, $\quad$ $\forall u=(u_1,\dots,u_n), \varphi=(\varphi_1,\dots,\varphi_n) \in X$.
		\end{center}

		Furthermore, $\Psi^\prime$ is compact from $X$ to $X^{\star}$.
	\end{lemma}

	\begin{proof}
	We start by showing that the functional $F$ is well defined on $X$. Indeed, for all $u=(u_1,\dots,u_n) \in X$, we have in virtue of $(\mathcal{F}_1)$ and $(\mathcal{F}_2)$,

		\begin{center}
			$\displaystyle \begin{aligned} F(x,u) & = \sum_{i=1}^n \int_{0}^{u_i} \frac{\partial F}{\partial s}(x, u_1, \dots, s, \dots, u_n) ds + F(x, 0, \dots, 0) dx \\ & \leq c_{1}\left[\sum_{i=1}^{n}\left(\sum_{j=1}^{n} b_{i j}(x)\left|u_{j}(x)\right|^{\mu_{i j}-1}\left|u_{i}(x)\right|\right)\right], \end{aligned}$
		\end{center}

		Then,

		\begin{center}
			$\displaystyle \int_{\mathbb{R}^{N}} F\left(x, u_{1}, \ldots, u_{n}\right) d x \leq c_{2}\left[\sum_{i=1}^{n}\left(\int_{\mathbb{R}^{N}} \sum_{j=1}^{n} b_{i j}(x)\left|u_{j}(x)\right|^{\mu_{i j}-1}\left|u_{i}(x)\right| d x\right)\right].$
		\end{center}

		If we consider the fact that $W^{1, \gamma(x)}\left(\mathbb{R}^{N}\right) \hookrightarrow L^{\mu(x)}\left(\mathbb{R}^{N}\right)$, for $\mu(x)>1$, then there exists $c>0$ such that

		\begin{center}
			$\displaystyle ||u|^{\mu}|_{\gamma(x)} = |u|_{\mu \gamma(x)}^{\mu} \leq c \| u \|_{\gamma(x)}^{\mu}$,
		\end{center}

		and if we apply Propositions \ref{prop1},\ref{prop2} and \ref{prop5} and take $b_{i i} \in L^{\alpha_{i}(x)}, b_{i j} \in L^{\alpha_{i j}(x)}$ if $i \neq j$, then we have

		\begin{align}
			\displaystyle  \int_{\mathbb{R}^{N}} F(x, u_{1}, \ldots, u_{n}) & \leq c_{3} \Big[\sum_{i=1}^{n} \Big(\sum_{j=1}^{n}|b_{i j}|_{\alpha_{i j}(x)}|| u_{j}|^{\mu_{i j}-1}|_{\gamma_{j}^{\partial}(x)}|u_{i}|_{\gamma_{i}^{\partial}(x)}\Big)\Big] \quad \quad (3.1) \\ & \leq c_{3} \Big[\sum_{i=1}^{n} \Big(\sum_{j=1}^{n}|b_{i j}|_{\alpha_{i j}(x)} |u_{j}|_{(\mu_{i j}-1) \gamma_{j}^{\partial}(x)}^{\mu_{i j}-1} |u_{i}|_{\gamma_{i}^{\partial}(x)}\Big)\Big] \quad  \quad (3.2) \\
				& \leq c_{3} \Big[\sum_{i=1}^{n} \Big(\sum_{j=1}^{n}|b_{i j}|_{\alpha_{i j}(x)}\|u_{j}\|_{\gamma_{j}(x)}^{\mu_{i j}-1}\|u_{i}\|_{\gamma_{i}(x)}\Big)\Big] < \infty.  \quad \label{(3.3)} \quad(3.3)
			\end{align}.
		\

		Hence, $\Psi$ is well defined. Moreover, one can easily see that $\Psi^{\prime}$ is also well defined on $X$. Indeed, using $(\mathcal{F}_2)$ for all $\varphi=\left(\varphi_{1}, \dots, \varphi_{n}\right) \in X$, we have

		\begin{center}
			$\displaystyle \begin{aligned} \Psi^{\prime}(u) \varphi &=\sum_{i=1}^{n} \int_{\mathbb{R}^{N}} \frac{\partial F}{\partial u_{i}}\left(x, u_{1}, \ldots, u_{n}\right) \varphi_{i} d x \\ & \leq \sum_{i=1}^{n}\left(\sum_{j=1}^{n} b_{i j}(x)\left|u_{j}(x)\right|^{\mu_{i j}-1}\right)\left|\varphi_{i}(x)\right| dx. \end{aligned}$
		\end{center}

		Following HÃ¶lder inequality, we obtain

		\begin{center}
			$\displaystyle \Psi^{\prime}(u) \varphi \leq c_{4} \Big[\sum_{i=1}^{n} \Big(\sum_{j=1}^{n}|b_{i j}|_{\alpha_{i j}(x)}|| u_{j}|^{\mu_{i j}-1}|_{\gamma_{j}^{\partial}(x)}|\varphi_{i}|_{\gamma_{i}^{\partial}(x)}\Big)\Big]$.
		\end{center}

		The above propositions yield

		\begin{center}
			$\displaystyle \Psi^{\prime}(u) \varphi \leq c\left[\sum_{i=1}^{n}\left(\sum_{j=1}^{n}\left|b_{i j}\right|_{\beta_{i j}(x)}\left\|u_{j}\right\|_{\gamma_{j}(x)}^{\mu_{i j}-1}\left\|\varphi_{i}\right\|_{\gamma_{i}(x)}\right)\right]<\infty$.
		\end{center}

		Now let us show that $\Psi$ is differentiable in the sense of Frechet, that is, for fixed $u=\left(u_{1}, \ldots, u_{n}\right) \in X$ and given $\varepsilon>0$, there must be a $\delta=\delta_{\varepsilon, u_{1}, \ldots, u_{n}}>0$ such that

		\begin{center}
			$\displaystyle \left|\Psi\left(u_{1}+\varphi_{1}, \ldots, u_{n}+\varphi_{n}\right)-\Psi\left(u_{1}, \ldots, u_{n}\right)-\Psi^{\prime}\left(u_{1}, \ldots, u_{n}\right)\left(\varphi_{1}, \ldots, \varphi_{n}\right)\right| \leq \varepsilon \sum_{i=1}^{n}\left(\left\|\varphi_{i}\right\|_{\gamma_{i}(x)}\right),$
		\end{center}

		for all $\varphi=\left(\varphi_{1}, \ldots, \varphi_{n}\right) \in X$ with $\displaystyle \sum_{i=1}^{n}\left(\left\|\varphi_{i}\right\|_{\gamma_{i}(x)}\right) \leq \delta$.

		Let $B_{R} = \Big\{ x \in \mathbb{R}^N: |x| < R \Big\}$ be the ball of radius $R$ which is centered at the origin of $\mathbb{R}^{N}$ and denote $B_{R}^{\prime}=\mathbb{R}^{N} \backslash B_{R}$. Moreover, let us define the functional $\Psi_{R}: \displaystyle \prod_{i=1}^{n} \Big(W_{w_{i}}^{1, p_{i}(x)}\left(B_{R}\right)\cap W_{w_{i}}^{1, \gamma_{i}(x)}\left(B_{R}\right)\Big) \longrightarrow \mathbb{R}$ as follows:

		\begin{center}
			$\displaystyle \Psi_{R}(u)=\int_{B_{R}} F\left(x, u_{1}(x), \ldots, u_{n}(x)\right) dx, \quad \forall u=(u_1,\dots,u_n)$.
		\end{center}

		If we consider $(\mathcal{F}_1)$ and $(\mathcal{F}_2)$, it is easy to see that $\displaystyle \Psi_{R} \in C^{1}\Big(\prod_{i=1}^{n}  W_{w_{i}}^{1, p_{i}(x)}\left(B_{R}\right)\cap W_{w_{i}}^{1, \gamma_{i}(x)}\left(B_{R}\right)\Big)$, and in addition for all $\displaystyle \varphi=\left(\varphi_{1}, \ldots, \varphi_{n}\right) \in \prod_{i=1}^{n}  \Big(W_{w_{i}}^{1, p_{i}(x)}\left(B_{R}\right)\cap W_{w_{i}}^{1, \gamma_{i}(x)}\left(B_{R}\right)\Big)$, we have

		\begin{center}
			$\displaystyle \Psi_{R}^{\prime}(u) \varphi=\sum_{i=1}^{n} \int_{B_{R}} \frac{\partial F}{\partial u_{i}}\left(x, u_{1}(x), \ldots, u_{n}(x)\right) \varphi_{i}(x) dx$.
		\end{center}

		Also as we know, the operator $\displaystyle \Psi_{R}^{\prime}: \prod_{i=1}^{n}  \Big(W_{w_{i}}^{1, p_{i}(x)}\left(B_{R}\right)\cap W_{w_{i}}^{1, \gamma_{i}(x)}\left(B_{R}\right)\Big) \longrightarrow \Big(\prod_{i=1}^{n}  \Big(W_{w_{i}}^{1, p_{i}(x)}\left(B_{R}\right)\cap W_{w_{i}}^{1, \gamma_{i}(x)}\left(B_{R}\right)\Big)\Big)^{\star}$ is compact \cite{Fan3}. Then, for all $u=\left(u_{1}, \ldots, u_{n}\right), \varphi=\left(\varphi_{1}, \ldots, \varphi_{n}\right) \in X$, we can write

		\begin{center}
			$\displaystyle \begin{aligned} & \left|\Psi\left(u_{1}+\varphi_{1}, \ldots, u_{n}+\varphi_{n}\right)-\Psi\left(u_{1}, \ldots, u_{n}\right)-\Psi^{\prime}\left(u_{1}, \ldots, u_{n}\right)\left(\varphi_{1}, \ldots, \varphi_{n}\right)\right| \\
		 & \leq\left|\Psi_{R}\left(u_{1}+\varphi_{1}, \ldots, u_{n}+\varphi_{n}\right)-\Psi_{R}\left(u_{1}, \ldots, u_{n}\right)-\Psi_{R}^{\prime}\left(u_{1}, \ldots, u_{n}\right)\left(\varphi_{1}, \ldots, \varphi_{n}\right)\right| \\
		&\quad + \Big| \int_{B_{R}^{\prime}}\left(F\left(x, u_{1}+\varphi_{1}, \ldots, u_{n}+\varphi_{n}\right)-F\left(x, u_{1}, \ldots, u_{n}\right)-\right. \left.\sum_{i=1}^{n} \int_{B_{R}} \frac{\partial F}{\partial u_{i}}\left(x, u_{1}, \ldots, u_{n}\right) \varphi_{i}\right) d x \Big| .
	\end{aligned}$
		\end{center}

		By virtue of mean-value theorem, there exist $\left.\xi_{1}, \ldots, \xi_{n} \in \right ]0,1[$ such that

		\begin{center}
			$\displaystyle \begin{gathered} \left| \int_{B_{R}^{\prime}}\left(F\left(x, u_{1}+\varphi_{1}, \ldots, u_{n}+\varphi_{n}\right)-F\left(x, u_{1}, \ldots, u_{n}\right)\right) - \sum_{i=1}^{n} \int_{B_{R}} \frac{\partial F}{\partial u_{i}}\left(x, u_{1}, \ldots, u_{n}\right) \varphi_{i} dx \right| \\ = \left|\int_{B_{R}^{\prime}}\left(\sum_{i=1}^{n} \frac{\partial F}{\partial u_{i}}\left(x, u_{1}, \ldots, u_{i}+\xi_{i} \varphi_{i}, \ldots, u_{n}\right) \varphi_{i}-\sum_{i=1}^{n} \frac{\partial F}{\partial u_{i}}\left(x, u_{1}, \ldots, u_{n}\right)\right) dx \right| .\end{gathered}$
		\end{center}

		Using the condition $(\mathcal{F}_2)$, we have

		\begin{center}
			$\displaystyle \left|\int_{B_{R}^{\prime}}\left(F\left(x, u_{1}+\varphi_{1}, \ldots, u_{n}+\varphi_{n}\right)-F\left(x, u_{1}, \ldots, u_{n}\right)\right)-\sum_{i=1}^{n} \int_{B_{R}} \frac{\partial F}{\partial u_{i}}\left(x, u_{1}, \ldots, u_{n}\right) \varphi_{i} d x\right| \leq\left|\sum_{i=1}^{n}\left(\sum_{j=1}^{n} \int_{B_{R}^{\prime}} b_{i j}(x)\left(\left|u_{j}+\xi_{j} \varphi_{j}\right|^{\mu_{i j}-1}-\left|u_{j}\right|^{\mu_{i j}-1}\right) \varphi_{i} d x\right)\right|.$
		\end{center}

		Using the elementary inequality $|a+b|^{s} \leq 2^{s-1}\left(|a|^{s}+|b|^{s}\right)$ for $a, b \in \mathbb{R}^{N}$, we can write

		\begin{center}
			$\displaystyle \begin{gathered} \leq \sum_{i=1}^{n}\left(\sum _ { j = 1 } ^ { n } \left(\left(2^{\mu_{i j}-1}-1\right) \int_{B_{R}^{\prime}} b_{i j}(x)\left|u_{j}\right|^{\mu_{i j}-1}\left|\varphi_{i}\right| dx \right.\left.+\left(\xi_{j} 2\right)^{\mu_{i j}-1} \int_{B_{R}^{\prime}} b_{i j}(x)\left|\varphi_{j}\right|^{\mu_{i j}-1}\left|\varphi_{i}\right| d x\right)\right) \end{gathered}$.
		\end{center}

		Then, applying Propositions \ref{prop2}, \ref{prop3} and \ref{prop5}, we have

		\begin{center}
			$\displaystyle \leq \sum_{i=1}^{n} c\left(\sum_{j=1}^{n}\left(\left|b_{i j}(x)\right|_{\alpha_{i j}}\left\|u_{j}\right\|_{\gamma_{1}^{\partial}(x)}^{\mu_{i j}-1}+\left|b_{i j}(x)\right|_{\alpha_{i j}}\left\|\varphi_{j}\right\|_{\gamma_{j}^{\partial}(x)}^{\mu_{i j}-1}\right)\right)\left\|\varphi_{i}\right\|_{\gamma_{i}(x)}$
		\end{center}

		and by the fact that

		\begin{center}
			$\displaystyle \left|b_{i i}(x)\right|_{L^{\alpha_{i}}\left(B_{R}^{\prime}\right)} \longrightarrow 0$,
		\end{center}

		and

		\begin{center}
			$\displaystyle \left|b_{i j}(x)\right|_{L^{\alpha_{i}}\left(B_{R}^{\prime}\right)} \longrightarrow 0$,
		\end{center}

		for all $1 \leq i, j \leq n$, as $R \rightarrow \infty$, and for $R$ sufficiently large, we obtain the estimate

		\begin{center}
			$\displaystyle \Big| \int_{B_{R}^{\prime}}\left(F\left(x, u_{1} + \varphi_{1}, \ldots, u_{n} + \varphi_{n}\right) - F\left(x, u_{1}, \ldots, u_{n}\right) - \sum_{i=1}^{n} \frac{\partial F}{\partial u_{i}}\left(x, u_{1}, \ldots, u_{n}\right) \varphi_{i}\right) dx \Big| \leq \varepsilon \sum_{i=1}^{n}\left(\left\|\varphi_{i}\right\|_{\gamma_{i}(x)}\right) $.
		\end{center}

	It only remains to show that $\Psi^{\prime}$ is continuous on $X$. Let $u_{m}=\left(u_{m,1}, \ldots, u_{m,n}\right)$ be such that $u_{m} \rightarrow u$ as $m \rightarrow \infty$. Then, for $\varphi = \left(\varphi_{1}, \ldots, \varphi_{n}\right) \in X$, we have

		\begin{center}
			$\displaystyle \begin{aligned} \left|\Psi^{\prime}\left(u_{m}\right) \varphi - \Psi^{\prime}(u) \varphi \right| & \leq \left|\Psi_{R}^{\prime}\left(u_{m}\right) \varphi -\Psi_{R}^{\prime}(u) \varphi \right| \\ & + \sum_{i=1}^{n} \int_{B_{R}^{\prime}}\left|\left(\frac{\partial F}{\partial u_{i}}\left(x, u_{m,1}, \ldots, u_{m,n}\right) \varphi_{i}-\frac{\partial F}{\partial u_{i}}\left(x, u_{1}, \ldots, u_{n}\right) \varphi_{i}\right) dx \right| \end{aligned}$
		\end{center}

		Since $\Psi_{R}^{\prime}$ is continuous on $\displaystyle \prod_{i=1}^{n} \Big(W_{w_{i}}^{1, p_{i}(x)}\left(B_{R}\right)\cap W_{w_{i}}^{1, \gamma_{i}(x)}\left(B_{R}\right)\Big)$ (see \cite{Fan3}), we have

		\begin{center}
			$\displaystyle \left|\Psi_{R}^{\prime}\left(u_{m}\right) \varphi -\Psi_{R}^{\prime}(u) \varphi\right| \longrightarrow 0\quad \text{ as } m \rightarrow \infty$.
		\end{center}

		Now, using $(\mathcal{F}_2)$ once again and taking into account that the other terms on the right-hand side of the above inequality tend to zero, we conclude that $\Psi^{\prime}$ is continuous on $X$. \\

		As for the compactness of $\Psi^{\prime}$, let $u_{m}=\left(u_{m,1}, \ldots, u_{m,n}\right)$ be a bounded sequence in $X$. Then, there exists a subsequence (we also denote it as $u_{m}=\left(u_{m,1}, \ldots, u_{m,n}\right)$) which converges weakly in $X$ to $u=\left(u_{1}, \ldots, u_{n}\right) \in X$. Then, if we use the same arguments as above, we have

		\begin{center}
			$\displaystyle \begin{aligned} \left|\Psi^{\prime}\left(u_{m}\right) \varphi -\Psi^{\prime}(u) \varphi\right| & \leq \left|\Psi_{R}^{\prime}\left(u_{m}\right) \varphi -\Psi_{R}^{\prime}(u) \varphi \right| \\ & + \sum_{i=1}^{n} \int_{B_{R}^{\prime}}\left|\left(\frac{\partial F}{\partial u_{i}}\left(x, u_{m,1}, \ldots, u_{m,n} \right) \varphi_{i}-\frac{\partial F}{\partial u_{i}}\left(x, u_{1}, \ldots, u_{n}\right) \varphi_{i}\right) d x\right|. \end{aligned}$
		\end{center}

		Since the restriction operator is continuous, we have $u_{m} \rightharpoonup u$ in  $\displaystyle \prod_{i=1}^{n} \Big(W_{w_{i}}^{1, p_{i}(x)}\left(B_{R}\right)\cap W_{w_{i}}^{1, \gamma_{i}(x)}\left(B_{R}\right)\Big)$. Because of the compactness of $\Psi^{\prime}$, the first expression on the right-hand side of the inequality tends to 0 , as $m \longrightarrow \infty$, and, as above, for sufficiently large $R$ we obtain

		\begin{center}
			$\displaystyle \sum_{i=1}^{n} \int_{B_{R}^{\prime}}\left|\left(\frac{\partial F}{\partial u_{i}}\left(x,u_{m,1}, \ldots, u_{m,n}\right) \varphi_{i}-\frac{\partial F}{\partial u_{i}}\left(x, u_{1}, \ldots, u_{n}\right) \varphi_{i}\right) \right| dx \longrightarrow 0.$
		\end{center}

		This implies $\Psi^{\prime}$ is compact from $X$ to $X^{\star}$.
	\end{proof}

	\begin{theorem}\label{Th3.3}
		Under assumptions $(\mathcal{F}_1)-(\mathcal{F}_3)$, the system \eqref{s1.1} admits at least three distinct weak solutions in $X$ for each $\lambda$ that belongs to the interval

		\begin{center}
			$\displaystyle \Lambda_{r}= \Bigg] \frac{\displaystyle m_1 \sum_{i=1}^n \max \left\lbrace \| z_i\|_{p_i(x),w_i}^{p_i^-} ,\| z_i\|_{p_i(x),w_i}^{p_i^+} \right\rbrace + \mathcal{H}(k_i^3) \max \left\lbrace \| z_i\|_{q_i(x),w_i}^{q_i^-} ,\| z_i\|_{q_i(x),w_i}^{q_i^+} \right\rbrace}{\displaystyle \int_{\mathbb{R}^N} F(x, z_1, \dots, z_n) dx} , \frac{r}{\displaystyle \int_{\mathbb{R}^N} \sup_{(\xi_1, \dots, \xi_n) \in K(\frac{sr}{m_0^{\star}})} F(x, \xi_1, \dots, \xi_n) dx} \Bigg[.$
		\end{center}
	\end{theorem}

	\begin{proof}
		By Lemma \ref{lem.3.1}, $\Phi$ is coercive, and by definitions of $\Phi$ and $\Psi$ and from hypothesis $(\mathcal{F}_1)$, we have $\Phi(0,\dots,0) = \Psi(0, \dots, 0) = 0$. Moreover, the required hypothesis $\Phi(\bar x) > r$ follows from $(C_1)$ and the definition of $\Phi$ by choosing $\bar x = (z_1, \dots, z_n)$.	On the contrary, by applying Proposition \ref{prop4} for $u=(u_1, \dots, u_n) \in X$, we have

			\begin{multline}\label{inq32}
				\displaystyle \frac{1}{s} \Bigg( \sum_{i=1}^n   \Big( \min\Big\{|u_i|_{p_i^{\partial}}^{p_i^-}, |u_i|_{p_i^{\partial}}^{p_i^+}\Big\} + \mathcal{H}(\kappa_i^3) \cdot \min\Big\{ |u_i|_{q_i^{\partial}}^{q_i^-}, |u_i|_{q_i^{\partial}}^{q_i^+} \Big\} \Big) \Bigg) \\ \leq \sum_{i=1}^n \frac{1}{p_i^+}  \Big(\min\Big\{\| u_i\|_{p_i(x),w_i}^{p_i^-} , \| u_i\|_{p_i(x),w_i}^{p_i^+}\Big\}  +\mathcal{H}(k_{i}^3) \cdot \min\Big\{\| u_i\|_{q_i(x),w_i}^{q_i^-}, \| u_i\|_{q_i(x),w_i}^{q_i^+}\Big\} \Big) .
				\end{multline}

		with $\displaystyle s = \max_{1 \leq i \leq n}\Big\{p_i^+ \displaystyle \min \{c_{p_i(x)}^{p_i^{-}}, c_{p_i(x)}^{p_i^{+}} \}, p_i^+ \displaystyle \min \{c_{q_i(x)}^{q_i^{-}}, c_{q_i(x)}^{q_i^{+}} \}\Big \} $, defined in $\eqref{defs}$. \\

		Now, from  \eqref{inq32}, we obtain for $r > 0$

		\begin{center}
			$\displaystyle \begin{aligned} \Phi^{-1}(]-\infty,r]) & = \Big\{u = (u_1, \dots, u_n) \in X: \Phi(u_1, \dots, u_n) \leq r \Big\} \\
				& \subseteq \Big\{u = (u_1, \dots, u_n) \in X: \displaystyle \sum_{i=1}^n  \frac{m_i^0\min\{\displaystyle \min\{k_{1i}^0,k_{2i}^0\}, \displaystyle \min \{k_{1i}^2,k_{2i}^2\}\}}{p_i^+\displaystyle \max\{\beta_{1i}, \beta_{2i}\}} \\ &
				\quad \cdot \Big(\min\{\| u_i\|_{p_i(x),w_i}^{p_i^-} , \| u_i\|_{p_i(x),w_i}^{p_i^+}\} + \mathcal{H}(k_{i}^3) \cdot \min\{\| u_i\|_{q_i(x),w_i}^{q_i^-}, \| u_i\|_{q_i(x),w_i}^{q_i^+}\} \Big) \leq r \Big\} \\
				& \subseteq \displaystyle\Big\{\frac{\displaystyle m_0 \min_{1 \leq i \leq n}\{\kappa_i\} }{\displaystyle\max_{1 \leq i \leq n}\Big\{p_i^+ \displaystyle \min \{c_{p_i(x)}^{p_i^{-}}, c_{p_i(x)}^{p_i^{+}} \}, p_i^+ \displaystyle \min \{c_{q_i(x)}^{q_i^{-}}, c_{q_i(x)}^{q_i^{+}} \}\Big \}} \sum_{i=1}^{n} \Big( \min\{ |\xi_i|_{p_i^{\partial}(x)}^{{p_i^{-}}}, |\xi_i|_{p_i^{\partial}(x)}^{{p_i^{+}}} \} \\
				& + \mathcal{H}(k_i^3) \min \{ |\xi_i|_{q_i^{\partial}(x)}^{{q_i^{-}}}, |\xi_i|_{q_i^{\partial}(x)}^{{q_i^{+}}} \} \Big) \leq r \Big\} \\
				 & = \displaystyle K\Big(\frac{sr}{m_0^{\star}}\Big) \end{aligned}$
		\end{center}

		where $K(t)$ is defined in $\eqref{defK}$, and $\kappa_i$ and $m_0^{\star}$ are defined as follows

		\begin{center}
			$\displaystyle \begin{aligned} \kappa_i & = \displaystyle \frac{\min\{\displaystyle \min\{k_{1i}^0,k_{2i}^0\}, \displaystyle \min \{k_{1i}^2,k_{2i}^2\}\}}{\displaystyle \max\{\beta_{1i}, \beta_{2i}\}},\\
				 m_0^{\star} & = \displaystyle m_0 \min_{1 \leq i \leq n}\{\kappa_i\} .\end{aligned}$
		\end{center}

Then, we get

		$\displaystyle \begin{aligned} \sup_{(u_1, \dots, u_n) \in \Phi^{-1}(]-\infty, r])} \Psi (u_1, \dots, u_n) & = \sup_{(\xi_1, \dots, \xi_n) \in K(\frac{sr}{m_0^{\star}})} \int_{\mathbb{R}^N} F(x, u_1, \dots, u_n) dx \\ & \leq \int_{\mathbb{R}^N} \sup_{(\xi_1, \dots, \xi_n) \in K(\frac{sr}{m_0^{\star}})} F(x, \xi_1, \dots, \xi_n) dx .\end{aligned}$ \\

		Therefore from $(C_2)$, we have

		$\displaystyle \begin{aligned} \sup_{(u_1, \dots, u_n) \in \Phi^{-1}(]-\infty, r])} \Psi (u_1, \dots, u_n) & \leq r \frac{\displaystyle \int_{\mathbb{R}^N} F(x, z_1, \dots, z_n) dx}{\displaystyle m_1 \sum_{i=1}^n \Big(\max \left\lbrace \| z_i\|_{p_i(x),w_i}^{p_i^-} ,\| z_i\|_{p_i(x),w_i}^{p_i^+} \right\rbrace + \mathcal{H}(k_i^3) \max \left\lbrace \| z_i\|_{q_i(x),w_i}^{q_i^-} ,\| z_i\|_{q_i(x),w_i}^{q_i^+} \right\rbrace\Big) }
			\\ & \leq r\frac{\Psi(z_1, \dots, z_n)}{\Phi(z_1, \dots, z_n)}, \end{aligned}$

		 from which condition $(a_1)$ of Lemma \ref{lem1.1} follows.

		 In the next step, we shall prove that for each $\lambda>0$, the energy functional $E_\lambda = \Phi - \lambda \Psi$ is coercive. For all $u=(u_1, \dots, u_n) \in X$, under the assumptions $(\textbf{H}_2)$ and $(\textbf{H}_3)$ and by using inequality \eqref{(3.3)}, we have

		\begin{center}
			$\displaystyle \begin{aligned} E_\lambda (u) & = \sum_{i=1}^n \widehat{M}_i\Big(\mathcal{B}_i(u_i)\Big) - \lambda \int_{\mathbb{R}^N} F(x, u_1, \dots, u_n) dx \\
				& \geq \sum_{i=1}^n \frac{m_0}{p_i^+} \int_{\mathbb R^N} \Big(\mathcal{A}_{1i}(|\nabla u_i|^{p_i(x)}) + w_i(x)\mathcal{A}_{2i}(|u_i|^{p_i(x)}) \Big) dx \\
				&\quad - \lambda \int_{\mathbb{R}^N} \Big( \sum_{i=1}^n \int_0^{u_i} \frac{\partial F}{\partial s} (x, u_1, \dots, s, \dots, u_n) ds + F(x, 0, \dots, 0) dx \Big) \\
				 & \geq \sum_{i=1}^n \frac{m_0}{p_i^+} \int_{\mathbb R^N} \Big( \frac{1}{\beta_{1i}} a_{1i}\big( |\nabla u_i|^{p_i(x)} \big) |\nabla u_i|^{p_i(x)} + \frac{1}{\beta_{2i}}w_i(x) a_{2i}\big( |u_i|^{p_i(x)} \big) |u_i|^{p_i(x)} \Big) dx \\
				  & \quad- \lambda c_{3} \Big[\sum_{i=1}^{n} \Big(\sum_{j=1}^{n}|b_{i j}|_{\alpha_{i j}(x)}\|u_{j}\|_{\gamma_{j}(x)}^{\mu_{i j}-1}\|u_{i}\|_{\gamma_{i}(x)}\Big)\Big] \\
				  & \geq \sum_{i=1}^n  \frac{m_0}{p_i^+\max\{\beta_{1i}, \beta_{2i}\}} \int_{\mathbb R^N} \Big(a_{1i}\big( |\nabla u_i|^{p_i(x)} \big) |\nabla u_i|^{p_i(x)} + w_i(x) a_{2i}\big(|u_i|^{p_i(x)} \big) |u_i|^{p_i(x)} \Big) dx \\
				  & \quad - \lambda c_{3} \Big[\sum_{i=1}^{n} \Big(\sum_{j=1}^{n}|b_{i j}|_{\alpha_{i j}(x)}\|u_{j}\|_{\gamma_{j}(x)}^{\mu_{i j}-1}\|u_{i}\|_{\gamma_{i}(x)}\Big)\Big] \\
				  & \geq \sum_{i=1}^n \frac{m_0\min\big\{ \min \{k_{1i}^0 , k_{2i}^0\} , \min \{k_{1i}^2 , k_{2i}^2\} \big\}}{p_i^+\max\{\beta_{1i}, \beta_{2i}\}}  \int_{\mathbb R^N}\Big[ \Big(| \nabla u_i|^{p_i(x)}+w_i(x) |u_i|^{p_i(x)} \Big)  \\
				  & \quad   + \mathcal{H}(k_{i}^3) \Big(| \nabla u_i|^{q_i(x)}  + w_i(x) |u_i|^{q_i(x)} \Big)\Big]dx  - \lambda c_{3} \Big[\sum_{i=1}^{n} \Big(\sum_{j=1}^{n}|b_{i j}|_{\alpha_{i j}(x)}\|u_{j}\|_{\gamma_{j}(x)}^{\mu_{i j}-1}\|u_{i}\|_{\gamma_{i}(x)}\Big)\Big] \\
				& \geq \sum_{i=1}^n \frac{m_0 \kappa_i}{p_i^+} \cdot \Big(\rho_{p_i(x)w_i}(u_i) + \rho_{q_i(x)w_i}(u_i)\Big) - \lambda c_{3} \Big[\sum_{i=1}^{n} \Big(\sum_{j=1}^{n}|b_{i j}|_{\alpha_{i j}(x)}\|u_{j}\|_{\gamma_{j}(x)}^{\mu_{i j}-1}\|u_{i}\|_{\gamma_{i}(x)}\Big)\Big] \end{aligned}$
		\end{center}

		By using Young's inequality and \eqref{ineq1}, we obtain

		\begin{center}
			$\displaystyle \begin{aligned} E_\lambda (u) & \geq \sum_{i=1}^n \frac{m_0 \kappa_i}{p_i^+}  \Big(\|u\|_{p_i(x),w_i}^{p_i^-} + \mathcal{H}(k_{i}^3) \|u\|_{q_i(x),w_i}^{q_i^-} \Big)  - \lambda c_{3} \left[ \sum_{i=1}^{n} \left( \sum_{j=1}^{n} \left|b_{i j}\right|_{\alpha_{i j}(x)} \left(\frac{\mu_{i j}-1}{\mu_{i j}}\left\|u_{j}\right\|_{\gamma_{j}(x)}^{\mu_{i j}}\right. \right. \right.\left.\left.\left. +\frac{1}{\mu_{i j}}\left\|u_{i}\right\|_{\gamma_{i}(x)}^{\mu_{i j}}\right)\right) \right] \\
				 &  \geq \sum_{i=1}^n \left[ \frac{m_0 \kappa_i}{p_i^+}  \Big(\|u\|_{p_i(x),w_i}^{p_i^-} + \mathcal{H}(k_{i}^3) \|u\|_{q_i(x),w_i}^{q_i^-} \Big)  - \lambda c_{4}\left(\sum_{j=1}^{n}\left(\left|b_{i j}\right|_{\alpha_{i j}(x)}\left\|u_{j}\right\|_{\gamma_{j}(x)}^{\mu_{i j}}\right.\right.\right. \left.\left.\left.+\left|b_{i j}\right|_{\alpha_{i j}(x)}\left\|u_{i}\right\|_{\gamma_{i}(x)}^{\mu_{i j}}\right)\right)\right]. \end{aligned} $
		\end{center}

		This shows that $(\Phi - \lambda \Psi)(u) \longrightarrow +\infty$ as $\|u\| \longrightarrow +\infty$, since we have $\displaystyle 1 < \mu_{ij} < \inf_{x \in \mathbb{R}^N} \gamma_i(x)$. Therefore, it is confirmed that $\Phi - \lambda \Psi$ is coercive on $X$, for every parameter $\lambda>0$ in particular for every $\displaystyle \lambda \in \Lambda_r = ]\frac{\Phi(\bar x)}{\Psi(\bar x)}, \frac{r}{\displaystyle \sup_{\Phi(\bar x) \leq r} \Psi(\bar x)}[$,
		and then the hypothesis $(a_2)$ of  Lemma \ref{lem1.1}  is also proved. Now, all the hypotheses of Lemma \ref{lem1.1} are satisfied. Note that the solutions of the equation $\Phi^\prime (u) - \lambda \Psi^\prime (u)$ are exactly the weak solutions of \eqref{s1.1}. Thus, for each

		\begin{center}
			$\displaystyle \lambda \in \Bigg] \frac{\displaystyle m_1 \sum_{i=1}^n\Big( \max \left\lbrace \| z_i\|_{p_i(x),w_i}^{p_i^-} ,\| z_i\|_{p_i(x),w_i}^{p_i^+} \right\rbrace + \mathcal{H}(k_i^3) \max \left\lbrace \| z_i\|_{q_i(x),w_i}^{q_i^-} ,\| z_i\|_{q_i(x),w_i}^{q_i^+} \right\rbrace\Big)}{\displaystyle \int_{\mathbb{R}^N} F(x, z_1, \dots, z_n) dx} , \frac{r}{\displaystyle \int_{\mathbb{R}^N} \sup_{(\xi_1, \dots, \xi_n) \in K(\frac{sr}{m_0^{\star}})} F(x, \xi_1, \dots, \xi_n) dx} \Bigg[$
		\end{center}

		system \eqref{s1.1} admits at least three distinct weak solutions in $X$.
	\end{proof}