Solution 1 (6674.361 seconds): Discrete optimization
Break the trajectory into $N$ timesteps of length $\Delta t$. The time after the $i$-th timestep will be $t_i=i\Delta t$, with the starting time $t_0=0$ and the final time $t_N$. Given an acceleration of $\vec a_i$ during the $i$-th timestep, the velocity at time $t_i$ will be $\vec v_i=\sum_{k=1}^{i}\vec a_k\Delta t$ (assuming $\vec v_0=\vec 0$). The displacement during each timestep will be the average velocity times $\Delta t$: $\vec x_i-\vec x_{i-1}=(\vec v_i+\vec v_{i-1})\Delta t/2$. This means the total displacement will be:
$$
\begin{align}
\vec x_i &= \sum_{j=1}^i\left[\left(\sum_{k=1}^{j}\vec a_k\Delta t\right)+\left(\sum_{k=1}^{j-1}\vec a_k\Delta t\right)\right]\Delta t/2 \\
&= \sum_{k=1}^i[(2i-2k+1)\vec a_k]\Delta t^2/2
\end{align}
$$
(Again assuming that $\vec x_0=\vec 0$.)
I assert without proof that an optimal trajectory will pass through the corner. This means that we can separate the trajectory into two parts or legs. Break the first leg into $M=\lfloor N/2\rfloor$ equal timesteps $\Delta t_1$, and the second leg into $N-M$ equal timesteps $\Delta t_2$. Since the initial and final velocities are zero, we can write the velocity at $t_M$ as either $\vec v_M = \sum_{k=1}^{M}\vec a_k \Delta t_1 = -\sum_{k=M+1}^{N}\vec a_k \Delta t_2$. The displacement of each leg is known and can be calculated as above. Note that in both the velocity and displacement equations, $\Delta t_i$ can be factored out. Finally, we restrict $|\vec a_i|\le a_\textrm{max}$. This gives us the following system of equations:
Minimize: $t_N=M\Delta t_1+(N-M)\Delta t_2$
Subject to:
$$|\vec a_i| \le a_\textrm{max},\ i\in[1,N]\\\\$$
$$\begin{align}
\left(\frac{\Delta \vec v_1}{\Delta t_1}\right) &= \sum_{k=1}^{M}\vec a_k \\
\left(\frac{\Delta \vec v_2}{\Delta t_2}\right) &= \sum_{k=M+1}^{N}\vec a_k \\
\left(\frac{\Delta \vec x_1}{\Delta t_1^2}\right) &= \sum_{k=1}^{M}\left(M-k+\frac{1}{2}\right)\vec a_k \\
\left(\frac{\Delta \vec x_2}{\Delta t_2^2}\right) &= \sum_{k=M+1}^{N}\left(M-k+\frac{1}{2}\right)\vec a_k \\
\end{align}$$
$$
\left(\frac{\Delta \vec v_1}{\Delta t_1}\right)\Delta t_1 = \left(\frac{\Delta \vec v_2}{\Delta t_2}\right)\Delta t_2
$$
$$\begin{align}
\left(\frac{\Delta \vec x_1}{\Delta t_1^2}\right)\Delta t_1^2 &= (10^6,0) \\
\left(\frac{\Delta \vec x_2}{\Delta t_2^2}\right)\Delta t_2^2 &= (0,10^8) \\
\end{align}
$$
(The factor of $\vec a_k$ in the sum for $\Delta\vec x_2$ comes from the fact that we are summing in reverse from $N$ to $M+1$, so the factor is $(N-M)-j+1/2$ where $j=N-k+1$. This is equal to ($N-M-N+k-1+1/2=-M+k-1/2$), but since we are calculating backwards in time we need to negate the result, coincidentally giving us $M-k+1/2$. This is only possible since we know $\vec v_N=\vec 0$.)
The reason for factoring out $\Delta t_i$ is that this system has only six quadratic constraints, plus six linear constraints and $N$ convex cone constraints. This simplification lets a solver like SCIP handle the problem efficiently. Here is an example specification of the problem and a solution found by SCIP. And here is the solution in the format the validator expects; the validator output is:
Time: 6674.36
Final position & velocity: 1000048.6673564073 99999992.4373888 0.0033702956625456526 0.009344372255398525
Note my own calculations give a final position of (999 999.988, 99 999 999.984), well within the 1 meter tolerance, and a final velocity of (9.8×10-7, -7.9×10-7), which is 'substantially' less than 1. The position at $t_M$ is (999 999.982, -9.1×10-12), so I am technically cutting the corner, but only by 18 mm.
Solution 2 (6674.291 seconds): Calculus of Variations
Following the formulation in D.G. Stechert, "On the Use of the Calculus of Variations in Trajectory Optimization Problems," 1963 (pdf), we can set up our equations of motion in terms of $x$ and $y$:
$$
\begin{align}
\phi_1 &= \dot x - u = 0 \\
\phi_2 &= \dot y - v = 0 \\
\phi_3 &= \dot u - a\cos\theta = 0 \\
\phi_4 &= \dot v - a\sin\theta = 0 \\
\end{align}
$$
(Where $u$ and $v$ are the $x$- and $y$-velocities, respectively). We then set up our Lagrange multipliers:
$$
F = \lambda_1 \phi_1 + \lambda_2 \phi_2 + \lambda_3 \phi_3 + \lambda_4 \phi_4
$$
Compute derivatives:
$$
\begin{align}
F_x &= 0 & F_{\dot x} &= \lambda_1 \\
F_y &= 0 & F_{\dot x} &= \lambda_2 \\
F_u &= -\lambda_1 & F_{\dot x} &= \lambda_3 \\
F_v &= -\lambda_2 & F_{\dot x} &= \lambda_4 \\
&& F_\theta &= \lambda_3 a\sin\theta - \lambda_4 a\cos\theta \\
\end{align}
$$
And set up the Euler-Lagrange equations:
$$
\begin{align}
\dot\lambda_1 &= 0 \\
\dot\lambda_2 &= 0 \\
\dot\lambda_3 &= -\lambda_1 \\
\dot\lambda_4 &= -\lambda_2 \\
\end{align} \\
\lambda_3 a\sin\theta - \lambda_4 a\cos\theta = 0
$$
From these we can see that:
$$
\begin{align}
\lambda_1 &= c_1 \\
\lambda_2 &= c_2 \\
\lambda_3 &= c_3 - c_1 t \\
\lambda_4 &= c_4 - c_2 t \\
\end{align} \\
\theta = \arctan\left(\frac{c_4 - c_2 t}{c_3 - c_1 t}\right)
$$
We can show that with the substitutions $\theta_0=\arctan\left(-\frac{c_1}{c_2}\right)$, $t_0 = \frac{c_1 c_3 + c_2 c_4}{c_1^2 + c_2^2}$, and $\tau=\frac{c_1 c_4-c_2 c_3}{c_1^2+c_2^2}$, we have:
$$
\theta = \theta_0 + \arctan\left(\frac{t-t_0}{\tau}\right)
$$
In fact, if we look at the plots of $\theta$ for each leg of the discrete solution, we see that this is indeed the case:
Now we do get a closed-form expression for $x$, $y$, $u$, and $v$ by integrating the equations of motion with this arctan control input, however they are complicated enough that finding the optimal parameters will require numerical solution techniques.
I find that the optimal parameters for the two legs are:
$$
\begin{align}
t_1 &= 498.925 & t_2 &= 6175.367 \\
\theta_{0,1} &= 0.543648 & \theta_{0,2} &= 3.238171 \\
t_{0,1} &= 256.232 & t_{0,2} &= 3018.254 \\
\tau_1 &= 113.423 & \tau_2 &= 29.827 \\
\end{align}
$$
And here are some plots of the resulting trajectory. Note the second image is scaled down significantly in the vertical, since the course is 100 times longer in y than in x.
