This is a start, I'll add to this answer later today, or someone else is welcome to add an answer with more quantitative information.
I'll look at a solar sail in heliocentric orbit used to raise the orbital distance from the Sun. For a given mirror size and mass, it means you need roughly a 45 degree tilt to maximize the bang for the buck (delta-v per unit time per unit mass).
Let the $\mathbf{y}$ axis point towards the Sun and the $\mathbf{x}$ axis point in the direction of orbital motion (so $\mathbf{\hat{y}}, \mathbf{\hat{x}}$ are proxies for $\mathbf{\hat{r}}, \mathbf{\hat{\theta}}$ respectively) We can add scale to the problem later, so right now let's define the dimensionless parabola as
$$y=\frac{1}{2}x^2-1$$
with the focus at $(0, 0)$. We can already guess that the part of the mirror
with the "best" thrust/mass will be somewhere around 45 degrees (at about $(1, 0)$) but let's see where it is.
The sail will have a mass per unit length, so we need $ds$:
$$ds = \sqrt{dx^2 + dy^2} $$
Also needed is the Force per unit length $\frac{dF}{ds}$:
$$ \frac{dF}{ds} = \frac{dF}{dx} \frac{dx}{ds} = \frac{dx \sin(2 \theta)}{\sqrt{dx^2 + dy^2}} $$
The plot highlights the 95% of maximum region:

lower peak upper
x 0.53 0.71 0.93
θ 28.0 35.3 42.8
Optically that looks like this:

But if you used one big one, your solar collector (PV array) would be far from the sail. It's much better to have many small ones so that you don't add as much mass supporting the collectors accurately so far away.

import numpy as np
import matplotlib.pyplot as plt
halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
degs, rads = 180/pi, pi/180
x = np.linspace(0, 1.2, 2001)
y = 0.5*x**2-1
dx = x[1:]-x[:-1]
dy = y[1:]-y[:-1]
ds = np.sqrt(dx**2 + dy**2)
slope = dy/dx
theta = np.arctan2(dy, dx)
dFdx = np.sin(2.*theta)
dFds1 = dFdx * (dx/ds)
dFds = (dx * np.sin(2.*theta))/np.sqrt(dx**2 + dy**2)
nmaxi = np.argmax(dFds)
dFdsmaxi = dFds[nmaxi]
xmaxi = x[nmaxi]
thmax = degs * theta[nmaxi]
n95L = np.argmax(dFds>0.95*dFdsmaxi)
n95R = np.argmax(dFds[nmaxi:] < 0.95*dFdsmaxi) + nmaxi
print x[n95L], x[nmaxi], x[n95R]
print degs*theta[n95L], degs*theta[nmaxi], degs*theta[n95R]
if True:
plt.figure()
plt.subplot(2, 1, 1)
plt.plot(x[:-1], dFds, '-k')
plt.plot(x[:-1][n95L:n95R+1], dFds[n95L:n95R+1], '-r', linewidth=6)
plt.xlabel('x', fontsize=16)
plt.ylabel('dF/ds', fontsize=16)
plt.subplot(2, 1, 2)
plt.plot(degs*theta, dFds, '-k')
plt.plot(degs*theta[n95L:n95R+1], dFds[n95L:n95R+1], '-r', linewidth=6)
plt.xlabel('theta (deg)', fontsize=16)
plt.ylabel('dF/ds', fontsize=16)
plt.show()
"""
lower peak upper
x 0.53 0.71 0.93
th 28.0 35.3 42.8
"""
xx = x[n95L:n95R+1]
yy = y[n95L:n95R+1]
xxL = np.array([xx[0], xx[0], 0])
yyL = np.array([1, yy[0], 0])
xxR = np.array([xx[-1], xx[-1], 0])
yyR = np.array([1, yy[-1], 0])
O = np.zeros(1)
if True:
plt.figure
plt.plot(x, y, '--k', linewidth=0.5)
dx = xx[-1] - xx[0]
dy = yy[-1] - yy[0]
for n in range(1):
plt.plot(xx + n*dx, yy + n*dy, '-b', linewidth=6)
plt.plot(xxL + n*dx, yyL + n*dy, '-b')
plt.plot(xxR + n*dx, yyR + n*dy, '-b')
plt.plot(O + n*dx, O + n*dy, 'ok')
plt.xlim(-0.1, 1.1)
plt.ylim(-1.1, 0.1)
plt.show()
if True:
plt.figure
#plt.plot(x, y )
dx = xx[-1] - xx[0]
dy = yy[-1] - yy[0]
for n in range(5):
plt.plot(xx + n*dx, yy + n*dy, '-b', linewidth=6)
plt.plot(xxL + n*dx, yyL + n*dy, '-b')
plt.plot(xxR + n*dx, yyR + n*dy, '-b')
plt.plot(O + n*dx, O + n*dy, 'ok')
plt.xlim(-0.1, 2.1)
plt.ylim(-1.1, 1.1)
plt.show()