We concern a family
{
(
u
ε
,
v
ε
)
}
ε
>
0
\{(u_{\varepsilon },v_{\varepsilon })\}_{\varepsilon > 0}
of solutions of the Lane-Emden system on a smooth bounded convex domain
Ω
\Omega
in
R
N
\mathbb {R}^N
\[
{
−
Δ
u
ε
=
v
ε
p
a
m
p
;
in
Ω
,
−
Δ
v
ε
=
u
ε
q
ε
a
m
p
;
in
Ω
,
u
ε
,
v
ε
>
0
a
m
p
;
in
Ω
,
u
ε
=
v
ε
=
0
a
m
p
;
on
∂
Ω
,
\begin {cases} -\Delta u_{\varepsilon } = v_{\varepsilon }^p & \text {in } \Omega , \\ -\Delta v_{\varepsilon } = u_{\varepsilon }^{q_{\varepsilon }} & \text {in } \Omega , \\ u_{\varepsilon },\, v_{\varepsilon } > 0 & \text {in } \Omega , \\ u_{\varepsilon } = v_{\varepsilon } =0 & \text {on } \partial \Omega , \end {cases}
\]
for
N
≥
4
N \ge 4
,
max
{
1
,
3
N
−
2
}
>
p
>
q
ε
\max \{1,\frac {3}{N-2}\} > p > q_{\varepsilon }
and small
\[
ε
≔
N
p
+
1
+
N
q
ε
+
1
−
(
N
−
2
)
>
0.
\varepsilon ≔\frac {N}{p+1} + \frac {N}{q_{\varepsilon }+1} - (N-2) > 0.
\]
This system appears as the extremal equation of the Sobolev embedding
W
2
,
(
p
+
1
)
/
p
(
Ω
)
↪
L
q
ε
+
1
(
Ω
)
W^{2,(p+1)/p}(\Omega ) \hookrightarrow L^{q_{\varepsilon }+1}(\Omega )
, and is also closely related to the Calderón-Zygmund estimate. Under the natural energy condition, we prove that the multiple bubbling phenomena may arise for the family
{
(
u
ε
,
v
ε
)
}
ε
>
0
\{(u_{\varepsilon },v_{\varepsilon })\}_{\varepsilon > 0}
, and establish a detailed qualitative and quantitative description. If
p
>
N
N
−
2
p > \frac {N}{N-2}
, the nonlinear structure of the system makes the interaction between bubbles so strong, so the determination process of the blow-up rates and locations is completely different from that of the classical Lane-Emden equation. If
p
≥
N
N
−
2
p \ge \frac {N}{N-2}
, the blow-up scenario is relatively close to that of the classical Lane-Emden equation, and only single-bubble solutions can exist. Even in the latter case, we have to devise a new method to cover all
p
p
near
N
N
−
2
\frac {N}{N-2}
. We also deduce a general existence theorem that holds on any smooth bounded domains.