Let $B$ be a standard Brownian motion in $\mathbb R$. Define the variables $$\begin{align*} X &= B_1, & Y &= \int_0^1B_s\mathrm ds, & Z&= \int_0^1B_s^2\mathrm ds. \end{align*}$$ It is known, see below, that $(X,Y,Z)$ admits a smooth density.

Is is true that the density of $(X,Y,Z)$ is a Schwartz function?

I think there should be a simple enough argument, and I'm making it more complicated than it actually is.

Let me share two observations about this triple. First of all, all of its components belong to low order Wiener chaos: specifically, $X$ and $Y$ belong to the first Wiener chaos, and $Z$ to the second. This makes it easy enough to compute the Malliavin matrix ($\langle DX,DX\rangle$ and its peers):
$$ \Gamma
= \begin{pmatrix}
1 & 1/2 & 2\alpha \\
1/2 & 1/3 & 2\alpha-\beta \\
2\alpha & 2\alpha-\beta & 4\gamma
\end{pmatrix} $$
with
$$ \begin{align*} \alpha&=\int_0^1tB_t\mathrm dt, & \beta&=\int_0^1t^2B_t\mathrm dt, & \gamma&=\int_0^1\Big(\int_t^1B_s\mathrm ds\Big)^2\mathrm dt. \end{align*} $$
Using this expression, we get $\det\Gamma=(\gamma-\alpha^2)/3 - (\alpha-\beta)^2$, and it would suffice to prove that $\mathbb P[\det\Gamma\leq\varepsilon]$ decreases faster than any power of $\varepsilon$, or that $1/\det\Gamma$ has moments of all order (the former implies the latter). This would follow from Proposition 2.1.5 in Nualart's *The Malliavin Calculus and Related Topics* (p. 103), as suggested by Fabrice Baudoin in the comments.

As for the second observation, the variable in question is the evaluation at time 1 of a hypoelliptic process started at $(0,0,0)$: if one replaces every 1 by a $t$ in the above definitions, we get a process $(X_t,Y_t,Z_t)$ solution to the hypoelliptic equation $$\begin{align*} \mathrm dX_t &= \mathrm dB_t, & \mathrm dY_t &= X_t\mathrm dt,& \mathrm dZ_t&= X_t^2\mathrm dt \end{align*}$$ such that the value at time 1 is the variable described above. This shows that the density of $(X_1,Y_1,Z_1)$ is smooth. One can find more about this approach (with slightly different conventions) in section IV-1.5 (page 118) of my PhD thesis, but here are a few key insights. We can find a (not so nice) candidate for the characteristic function of $(X_t,Y_t,Z_t)$ (hence for $(X,Y,Z)$ as well) in the form $$\text{constant}\times\exp\big(\text{quadratic form in }(\xi_x,\xi_y)\big)$$ for $t$ and $\xi_z$ fixed. Specifically, $$ \mathbb E\big[\exp\big(-i(\xi_xX_t+\xi_yY_t+\xi_zZ_t)\big)\big] = \frac1{\sqrt{\mathrm{cah}}}\cdot \exp\left( - \frac t2\mathrm{tah}\cdot\xi_x^2 + \frac i{2\xi_z}\Big(1-\frac1{\mathrm{cah}}\Big)\cdot\xi_x\xi_y + \frac{it}{4\xi_z}(1-\mathrm{tah})\cdot\xi_y^2\right), $$ for $\mathrm{cah} = \cosh(\zeta)$, $\mathrm{sah} = \sinh(\zeta)/\zeta$, $\zeta = \sqrt{2i\xi_zt^2}$ (these are quantities of the form $f(\sqrt{(\cdots)})$ for $f$ holomorphic and even, so they are actually holomorphic), $\mathrm{tah} = \mathrm{cah}/\mathrm{sah}$. One can show that the zeros of $\alpha\mapsto\cosh(\sqrt\alpha)$ are real negative, so everything can be defined holomorphically in the above formula.

This is a nice candidate because it is *one* solution of the Feynman-Kac formula seen in Fourier space, in a rather strong sense. If it were the actual characteristic function, then the question would be equivalent to show that the candidate is itself a Schwartz function (in the dual space variables). However, this possibly equivalent result does not seem obvious to me either, since one would have to show that the real part of the quadratic form is bounded below, or something of the sort.