[proofplan]
**This proof is a sketch.** The proof has two directions. The forward direction $f \in H^1 \Rightarrow R_j f \in L^1$ uses the [atomic decomposition](/theorems/3177) and the standard two-region atomic estimate; this part is essentially complete. The reverse direction is sketch-grade and **uses two non-trivial cited results as black boxes**:
(B1) **Subharmonicity of $|F|^q$ for $q \ge (n-1)/n$** — for the harmonic vector $F = (u, v_1, \ldots, v_n)$ on $\mathbb{R}^{n+1}_+$ satisfying the generalised Cauchy–Riemann equations, the function $|F|^q$ is subharmonic for any $q \ge (n-1)/n$. This is a deep classical theorem due to Stein and Weiss; we cite [Stein–Weiss, *On the Theory of Harmonic Functions of Several Variables, II: Behaviour Near the Boundary*, Acta Math. 106 (1961), §3–4] for its proof, and **do not reprove it here**.
(B2) **Existence of the $L^{1/q}$ boundary trace of $|F|^q$** — under the hypotheses $u(\cdot, t), v_j(\cdot, t) \in L^1$ uniformly in $t$, the trace $|F(\cdot, 0)| := \lim_{t \to 0^+}|F(\cdot, t)|$ exists $\mathcal{L}^n$-a.e. and in $L^1$, and consequently $|F(\cdot, 0)|^q \in L^{1/q}(\mathbb{R}^n)$. This is also classical; see [Stein, *Harmonic Analysis: Real-Variable Methods, Orthogonality, and Oscillatory Integrals*, Ch. III] for a self-contained treatment.
Granted (B1) and (B2), the strategy is: starting from $f, R_j f \in L^1$, form the Poisson harmonic extensions $u, v_1, \ldots, v_n$ on $\mathbb{R}^{n+1}_+$ and observe that $F := (u, v_1, \ldots, v_n)$ satisfies the generalised Cauchy–Riemann system. Pick $q \in ((n-1)/n, 1)$. By (B1), $|F|^q$ is subharmonic, hence dominated by its Poisson extension; by (B2) the boundary trace exists in $L^{1/q}$; the Hardy–Littlewood maximal theorem on $L^{1/q}$ (which applies since $1/q > 1$) and the pointwise domination of Poisson by Hardy–Littlewood maximal then give $\sup_t |F(\cdot, t)| \in L^1$, hence $u^* := \sup_t |u(\cdot, t)| \in L^1$, hence $f \in H^1$. Full details of the cited inputs are in [Fefferman–Stein, *Acta Math.* 129 (1972)] and [Stein, *Harmonic Analysis*, Ch. III].
[/proofplan]
[step:Set up Riesz transforms and the harmonic extension]
The Riesz transforms are defined for $f \in \mathcal{S}(\mathbb{R}^n)$ by the Fourier multiplier
\begin{align*}
\widehat{R_j f}(\xi) = -i\, \frac{\xi_j}{|\xi|}\, \hat f(\xi), \quad j = 1, \ldots, n,
\end{align*}
or equivalently as the principal-value singular integral
\begin{align*}
R_j &: \mathcal{S}(\mathbb{R}^n) \to \mathcal{S}'(\mathbb{R}^n), \\
R_j f(x) &= c_n\, \mathrm{p.v.}\int_{\mathbb{R}^n} \frac{x_j - y_j}{|x - y|^{n+1}}\, f(y)\, d\mathcal{L}^n(y), \quad c_n := \Gamma((n+1)/2)/\pi^{(n+1)/2}.
\end{align*}
By the [$L^p$ boundedness of Riesz transforms](/theorems/???), $R_j: L^p(\mathbb{R}^n) \to L^p(\mathbb{R}^n)$ for $1 < p < \infty$, and $R_j: L^1(\mathbb{R}^n) \to L^{1,\infty}(\mathbb{R}^n)$. The Hardy space is
\begin{align*}
H^1(\mathbb{R}^n) := \{ f \in \mathcal{S}'(\mathbb{R}^n) : M_\varphi f \in L^1(\mathbb{R}^n) \},
\end{align*}
with $\varphi \in \mathcal{S}(\mathbb{R}^n)$ any fixed kernel with $\int\varphi\, d\mathcal{L}^n = 1$ ([independent of $\varphi$ by theorem 3173](/theorems/3173)). For the Poisson kernel
\begin{align*}
P_t(x) := c_n \frac{t}{(t^2 + |x|^2)^{(n+1)/2}}, \quad c_n = \Gamma((n+1)/2)/\pi^{(n+1)/2},
\end{align*}
the harmonic extension of $f$ to the upper half-space $\mathbb{R}^{n+1}_+ := \{(x,t) : x \in \mathbb{R}^n, t > 0\}$ is
\begin{align*}
u &: \mathbb{R}^{n+1}_+ \to \mathbb{R}, & u(x,t) &= (P_t * f)(x).
\end{align*}
For $f \in L^p(\mathbb{R}^n)$, $1 \le p < \infty$, $u$ is harmonic in $\mathbb{R}^{n+1}_+$, $u(\cdot, t) \to f$ in $L^p$ as $t \to 0^+$, and $\sup_{t>0}\|u(\cdot,t)\|_{L^p} = \|f\|_{L^p}$.
[/step]
[step:Forward direction — Riesz transforms map $H^1$ to $L^1$ via the atomic decomposition]
Suppose $f \in H^1(\mathbb{R}^n)$. By the [Coifman-Latter atomic decomposition](/theorems/3177), $f$ admits a decomposition
\begin{align*}
f = \sum_{k=1}^\infty \lambda_k\, a_k \quad \text{(convergence in $H^1$)},
\end{align*}
where each $a_k$ is an [$H^1$-atom](/theorems/???) supported in a ball $B_k$ with $\|a_k\|_\infty \le \mathcal{L}^n(B_k)^{-1}$ and $\int a_k\, d\mathcal{L}^n = 0$, and $\sum_k |\lambda_k| \asymp \|f\|_{H^1}$.
For each atom $a_k$ and each $j$, we claim
\begin{align*}
\|R_j a_k\|_{L^1} \le C_n,
\end{align*}
with $C_n$ depending only on $n$. The proof is the standard two-region split: let $B^* := 2 B_k$ denote the ball with double radius and the same centre as $B_k$. On $B^*$, by Cauchy-Schwarz and the $L^2$ boundedness of $R_j$,
\begin{align*}
\int_{B^*} |R_j a_k|\, d\mathcal{L}^n \le \mathcal{L}^n(B^*)^{1/2}\, \|R_j a_k\|_{L^2} \le C_n\, \mathcal{L}^n(B_k)^{1/2}\, \|a_k\|_{L^2} \le C_n\, \mathcal{L}^n(B_k)^{1/2}\, \mathcal{L}^n(B_k)^{1/2}\, \mathcal{L}^n(B_k)^{-1} = C_n.
\end{align*}
On $\mathbb{R}^n \setminus B^*$, use the cancellation $\int a_k\, d\mathcal{L}^n = 0$: writing $K_j(x) := c_n\, x_j/|x|^{n+1}$ for the Riesz kernel,
\begin{align*}
R_j a_k(x) = \int_{B_k} \bigl(K_j(x - y) - K_j(x - x_k)\bigr)\, a_k(y)\, d\mathcal{L}^n(y),
\end{align*}
where $x_k$ is the centre of $B_k$. The kernel difference is bounded for $x \notin 2B_k$ and $y \in B_k$ by $|y - x_k|^{1}\, |x - x_k|^{-(n+1)} \le r_k\, |x - x_k|^{-(n+1)}$, where $r_k = \operatorname{rad}(B_k)$, by the regularity estimate for the Riesz kernel. Therefore
\begin{align*}
|R_j a_k(x)| \le \int_{B_k} r_k\, |x - x_k|^{-(n+1)}\, |a_k(y)|\, d\mathcal{L}^n(y) \le r_k\, |x - x_k|^{-(n+1)}\cdot \mathcal{L}^n(B_k)\cdot \mathcal{L}^n(B_k)^{-1} = r_k\, |x - x_k|^{-(n+1)}.
\end{align*}
Integrating over $\{|x - x_k| > 2 r_k\}$ in polar coordinates (centred at $x_k$, with $s = |x - x_k|$ and surface element $\omega_{n-1}\,s^{n-1}\,d\mathcal{L}^1(s)$):
\begin{align*}
\int_{|x-x_k| > 2 r_k} r_k\,|x - x_k|^{-(n+1)}\,d\mathcal{L}^n(x) = r_k\,\omega_{n-1}\int_{2 r_k}^\infty s^{-(n+1)}\,s^{n-1}\,d\mathcal{L}^1(s) = r_k\,\omega_{n-1}\int_{2 r_k}^\infty s^{-2}\,d\mathcal{L}^1(s).
\end{align*}
The radial integral evaluates as $\int_{2 r_k}^\infty s^{-2}\,d\mathcal{L}^1(s) = [-s^{-1}]_{2 r_k}^\infty = (2 r_k)^{-1}$, so
\begin{align*}
\int_{|x-x_k| > 2 r_k} r_k\,|x - x_k|^{-(n+1)}\,d\mathcal{L}^n(x) = r_k\,\omega_{n-1}\cdot \frac{1}{2 r_k} = \frac{\omega_{n-1}}{2} =: C_n,
\end{align*}
where $\omega_{n-1}$ is the surface measure of $S^{n-1}$. Combining the two regions, $\|R_j a_k\|_{L^1} \le C_n$.
By the triangle inequality on the atomic decomposition,
\begin{align*}
\|R_j f\|_{L^1} \le \sum_k |\lambda_k|\, \|R_j a_k\|_{L^1} \le C_n\, \sum_k |\lambda_k| \le C'_n\, \|f\|_{H^1}.
\end{align*}
Summing over $j = 1, \ldots, n$ and combining with $\|f\|_{L^1} \le \|f\|_{H^1}$ ([since $H^1 \hookrightarrow L^1$](/theorems/???)) gives
\begin{align*}
\|f\|_{L^1} + \sum_{j=1}^n \|R_j f\|_{L^1} \le C(n)\, \|f\|_{H^1}.
\end{align*}
[/step]
[step:Reverse direction — set up the conjugate harmonic vector]
Now suppose $f \in L^1(\mathbb{R}^n)$ and $R_j f \in L^1(\mathbb{R}^n)$ for all $j = 1, \ldots, n$. Define the harmonic extensions
\begin{align*}
u, v_1, \ldots, v_n &: \mathbb{R}^{n+1}_+ \to \mathbb{R}, \\
u(x,t) &:= (P_t * f)(x), \\
v_j(x,t) &:= (P_t * R_j f)(x), \quad j = 1, \ldots, n.
\end{align*}
Each function is harmonic in $\mathbb{R}^{n+1}_+$ (the Poisson kernel $P_t$ produces harmonic extensions of its $L^1$ data). The vector $F := (u, v_1, \ldots, v_n): \mathbb{R}^{n+1}_+ \to \mathbb{R}^{n+1}$ satisfies the [generalised Cauchy-Riemann equations](/theorems/???):
\begin{align*}
\partial_t u + \sum_{j=1}^n \partial_{x_j} v_j &= 0, \\
\partial_{x_j} v_k - \partial_{x_k} v_j &= 0, \quad 1 \le j, k \le n, \\
\partial_t v_j - \partial_{x_j} u &= 0, \quad j = 1, \ldots, n.
\end{align*}
This is verified on the Fourier side: $\widehat{u(\cdot,t)}(\xi) = e^{-t|\xi|}\hat f(\xi)$ and $\widehat{v_j(\cdot,t)}(\xi) = -i(\xi_j/|\xi|)\, e^{-t|\xi|}\hat f(\xi)$, so
\begin{align*}
\partial_t u(\cdot, t)^\wedge(\xi) &= -|\xi|\, e^{-t|\xi|}\hat f(\xi), \\
\partial_{x_j} v_j(\cdot, t)^\wedge(\xi) &= i\xi_j\cdot (-i\xi_j/|\xi|)\, e^{-t|\xi|}\hat f(\xi) = (\xi_j^2/|\xi|)\, e^{-t|\xi|}\hat f(\xi),
\end{align*}
and summing over $j$ gives $\sum_j \xi_j^2/|\xi| = |\xi|$, so $\partial_t u + \sum_j \partial_{x_j} v_j = 0$. The other equations are verified similarly. Boundary values: $u(\cdot, t) \to f$ and $v_j(\cdot, t) \to R_j f$ in $L^1$ as $t \to 0^+$, and uniformly in $t$,
\begin{align*}
\sup_{t > 0}\|u(\cdot, t)\|_{L^1} \le \|f\|_{L^1}, \qquad \sup_{t > 0}\|v_j(\cdot, t)\|_{L^1} \le \|R_j f\|_{L^1},
\end{align*}
because $P_t \ge 0$ and $\|P_t\|_{L^1} = 1$.
[/step]
[step:Cite the subharmonic exponent inequality and the $L^{1/q}$ trace existence]
**Two black-box inputs are used here, both cited rather than reproved.**
(B1) **Stein–Weiss subharmonicity.** For the Cauchy–Riemann system of Step 3, [Stein–Weiss, *Acta Math.* 106 (1961), §3–4](/theorems/???) establishes that the function
\begin{align*}
G &: \mathbb{R}^{n+1}_+ \to [0, \infty), \\
G(x, t) &:= |F(x,t)|^q = \Bigl(u(x,t)^2 + \sum_{j=1}^n v_j(x,t)^2\Bigr)^{q/2}
\end{align*}
is **subharmonic on $\mathbb{R}^{n+1}_+$** for every $q \ge q_0 := (n-1)/n$. The exponent $q_0 < 1$ is sharp for the validity of the subharmonicity in this generality, and $q_0 < 1$ is precisely what allows the $H^1$ argument to close in the next step (since the Hardy–Littlewood maximal theorem on $L^{1/q}$ applies only when $1/q > 1$, i.e. $q < 1$). **We use this as a cited result and do not reprove it.**
(B2) **Existence of the $L^{1/q}$ boundary trace.** From $u(\cdot, t), v_j(\cdot, t) \in L^1$ uniformly in $t$ (Step 3), the trace
\begin{align*}
|F(\cdot, 0)| := \lim_{t \to 0^+}|F(\cdot, t)|
\end{align*}
exists $\mathcal{L}^n$-a.e. and in $L^1(\mathbb{R}^n)$, with $\||F(\cdot, 0)|\|_{L^1} \le \|f\|_{L^1} + \sum_j \|R_j f\|_{L^1}$ (since $u(\cdot, 0) = f$ and $v_j(\cdot, 0) = R_j f$ a.e. as boundary values of the Poisson extensions). Consequently,
\begin{align*}
|F(\cdot, 0)|^q \in L^{1/q}(\mathbb{R}^n), \quad \||F(\cdot, 0)|^q\|_{L^{1/q}}^{1/q} = \||F(\cdot, 0)|\|_{L^1}.
\end{align*}
**We cite this trace existence as a black box** ([Stein, *Harmonic Analysis*, Ch. III]).
**Consequence used below.** Granted (B1) and (B2), for any $q \in (q_0, 1)$, since $G \ge 0$ is subharmonic on $\mathbb{R}^{n+1}_+$ with $L^{1/q}$ boundary trace, the [Poisson representation for non-negative subharmonic functions](/theorems/???) yields
\begin{align*}
G(x, t) \le \bigl(P_t * G(\cdot, 0)\bigr)(x), \quad \text{i.e.} \quad |F(\cdot, t)|^q \le P_t * |F(\cdot, 0)|^q.
\end{align*}
[/step]
[step:Bound the Poisson maximal function of $f$ by $L^{1/q}$ data]
Pick $q \in (\tfrac{n-1}{n}, 1)$ — explicitly, $q := \tfrac{n-1}{n} + \tfrac{1}{2n} \in (q_0, 1)$ — so $1/q > 1$. The Poisson maximal function of $f$ satisfies
\begin{align*}
u^*(x) := \sup_{t > 0}|(P_t * f)(x)| = \sup_{t > 0}|u(x, t)| \le \sup_{t > 0}|F(x, t)|.
\end{align*}
Since $|F(x, t)| = (|F(x, t)|^q)^{1/q}$ and $a \mapsto a^{1/q}$ is monotone increasing on $[0, \infty)$,
\begin{align*}
\sup_{t > 0}|F(x, t)| = \Big(\sup_{t > 0}|F(x, t)|^q\Big)^{1/q}.
\end{align*}
\medskip
\noindent \textbf{Pointwise majorisation by the Hardy–Littlewood maximal function.} By Step 4,
\begin{align*}
|F(x, t)|^q \le (P_t * |F(\cdot, 0)|^q)(x).
\end{align*}
The Poisson kernel is dominated pointwise by a constant multiple of the Hardy–Littlewood centred maximal function: for any non-negative $g \in L^1_{\mathrm{loc}}(\mathbb{R}^n)$,
\begin{align*}
(P_t * g)(x) \le C(n)\,M g(x) \quad \text{for every } t > 0, \text{ for every } x \in \mathbb{R}^n,
\end{align*}
where $C(n)$ is a dimensional constant arising from the $L^1$-norm of the layer-cake decomposition of $P_t$ (this is the standard [pointwise domination of Poisson by Hardy–Littlewood maximal](/theorems/???)). Applying this with $g = |F(\cdot, 0)|^q$ and taking the supremum over $t > 0$:
\begin{align*}
\sup_{t > 0}|F(x, t)|^q \le C(n)\,M(|F(\cdot, 0)|^q)(x).
\end{align*}
\medskip
\noindent \textbf{$L^{1/q}$ boundedness of the Hardy–Littlewood maximal function.} Since $1/q > 1$, by the [Hardy–Littlewood maximal theorem](/theorems/???) the maximal operator $M$ is bounded on $L^{1/q}(\mathbb{R}^n)$ with operator norm $A_q = A(n, q)$:
\begin{align*}
\|M g\|_{L^{1/q}} \le A_q\,\|g\|_{L^{1/q}} \quad \text{for } g \in L^{1/q}(\mathbb{R}^n).
\end{align*}
\medskip
\noindent \textbf{The boundary data is in $L^{1/q}$.} The trace exists in $L^1$ (Step 4(B2)) with
\begin{align*}
\||F(\cdot, 0)|\|_{L^1} \le \|f\|_{L^1} + \sum_{j=1}^n \|R_j f\|_{L^1} =: D_f.
\end{align*}
For non-negative $G$, the identity $\|G\|_{L^{1/q}}^{1/q} = \big(\int G^{1/q}\,d\mathcal{L}^n\big)^q = \|G^{1/q}\|_{L^1}^q$ relates the two norms. Apply with $G = |F(\cdot, 0)|^q$ (so $G^{1/q} = |F(\cdot, 0)|$):
\begin{align*}
\big\||F(\cdot, 0)|^q\big\|_{L^{1/q}} = \big(\||F(\cdot, 0)|\|_{L^1}\big)^q = D_f^q.
\end{align*}
Hence $|F(\cdot, 0)|^q \in L^{1/q}(\mathbb{R}^n)$.
\medskip
\noindent \textbf{Closing the chain.} Combining the three steps above:
\begin{align*}
\big\|\sup_{t > 0}|F(\cdot, t)|^q\big\|_{L^{1/q}} \le C(n)\,\big\|M(|F(\cdot, 0)|^q)\big\|_{L^{1/q}} \le C(n)\,A_q\,\big\||F(\cdot, 0)|^q\big\|_{L^{1/q}} = C(n)\,A_q\,D_f^q.
\end{align*}
Apply the same identity $\|G\|_{L^{1/q}}^{1/q} = \|G^{1/q}\|_{L^1}^q$ in the *reverse direction*, this time with $G := \sup_{t > 0}|F(\cdot, t)|^q$ (so $G^{1/q} = \sup_{t > 0}|F(\cdot, t)|$ since $a \mapsto a^{1/q}$ is monotone), giving
\begin{align*}
\big\|\sup_{t > 0}|F(\cdot, t)|\big\|_{L^1}^q = \big\|\sup_{t > 0}|F(\cdot, t)|^q\big\|_{L^{1/q}} \le C(n)\,A_q\,D_f^q.
\end{align*}
Taking the $1/q$-th power:
\begin{align*}
\big\|\sup_{t > 0}|F(\cdot, t)|\big\|_{L^1} \le \big(C(n)\,A_q\big)^{1/q}\,D_f.
\end{align*}
Therefore
\begin{align*}
\|u^*\|_{L^1} \le \big\|\sup_{t > 0}|F(\cdot, t)|\big\|_{L^1} \le C(n, q)\,\Big(\|f\|_{L^1} + \sum_{j=1}^n \|R_j f\|_{L^1}\Big),
\end{align*}
where $C(n, q) := \big(C(n)\,A_q\big)^{1/q}$.
[/step]
[step:Combine the two directions]
Steps 2 (forward) and 5 (reverse, via Steps 3-4) together establish the equivalence of conditions: $f \in H^1(\mathbb{R}^n)$ iff $f \in L^1$ and all $R_j f \in L^1$. The constants from each direction are explicit in dimension $n$ (and the chosen exponent $q$ for the reverse direction, which can be fixed at $q := (n-1)/n + \tfrac{1}{2n} \in (q_0, 1)$). Combining,
\begin{align*}
\|f\|_{H^1(\mathbb{R}^n)} \asymp_n \|f\|_{L^1(\mathbb{R}^n)} + \sum_{j=1}^n \|R_j f\|_{L^1(\mathbb{R}^n)}.
\end{align*}
The Riesz transforms, originally defined on $L^1$ as bounded operators into $L^{1,\infty}$, extend continuously to $H^1 \to L^1$ by the forward direction; this extension agrees with the original on the dense subspace $L^1 \cap H^1$ because $L^{1,\infty}$ convergence and $L^1$ convergence agree on a sequence whose $L^1$ limit exists, completing the description of the operator domain.
[/step]